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1  Introduction 


Modern  space  surveillance  requires  fast  and  accurate  orbit  predictions  for  myriads  of  objects 
in  a  broad  range  of  Earth  orbits.  Conventional  Special  Perturbations  orbit  propagators, 
based  on  numerical  integration  of  the  osculating  equations  of  motion,  are  accurate  but  ex¬ 
tremely  slow  (typically  requiring  100  or  more  steps  per  satellite  revolution  to  give  good 
predictions).  Conventional  General  Perturbations  orbit  propagators,  based  on  fully  analyt¬ 
ical  orbit  theories  like  those  of  Brouwer,  are  faster  but  contain  large  errors  due  to  inherent 
approximations  in  the  theories.  New  orbit  generators  based  on  Semianalytic  Satellite  Theory 
(SST)  have  been  able  to  approach  the  accuracy  of  Special  Perturbations  propagators  and 
the  speed  of  General  Perturbations  propagators. 

SST  has  been  originated  by  P.  J.  Cefola  and  his  colleagues,  whose  names  are  in  the  refer¬ 
ences  at  the  end  of  this  document.  The  theory  is  scattered  throughout  the  listed  conference 
preprints,  published  papers,  technical  reports,  and  private  communications.  Our  purpose  in 
this  document  is  to  simplify,  assemble,  unify,  and  extend  the  theory. 

SST  represents  the  orbital  state  of  a  satellite  with  an  equinoctial  element  set  (aj,  •  •  •  a6). 
The  first  five  elements  aj,  •  •  •  a5  are  slowly  varying  in  time.  The  sixth  element  a6  is  the  mean 
longitude  A  and  so  is  rapidly  varying. 

SST  decomposes  the  osculating  elements  a,  into  mean  elements  a,  plus  a  small  remainder 
which  is  2 n  periodic  in  the  fast  variable: 

a,  =a,  +  j/t(a,,---a6,  t)  (1) 

(Here  we  use  hats  to  distinguish  the  elements  of  the  osculating  ellipse  from  the  elements  of 
the  averaging  procedure.  The  values  of  a  free  index  are  assumed  to  be  obvious  from  the 
context;  e.  g.  ,  here  i  can  have  the  values  1,  2,  3,  4,  5,  or  6,  so  (1)  represents  6  equations.) 
The  mean  elements  a,  are  governed  by  ordinary  differential  equations  of  the  form 

dat 

—  =  n6t6  +  Ai{au---a5,t)  (2) 

Here  t  is  the  time,  n  is  the  (mean)  mean  motion,  and  6#  is  the  Kronecker  delta  (i.  e.  , 
^16  =  ^26  =  ^36  =  a4v,  =  a56  =  0,  <$66  =  1).  The  short-periodic  variations  r/,  are  expressable  in 
Fourier  series  of  the  form 

OO 

Vi  =  £[<?/(«„ -«5,  t) cosjA  -f  5/(ctj, •  •  •  a5,  t) sinj’A]  (3) 

3-1 

Having  formulas  for  the  mean  element  rates  /!,,  we  can  integrate  the  mean  equations  (2) 
numerically  using  large  step  sizes  (typically  1  day  in  length).  The  formulas  for  the  Fourier 
coefficients  Cf  and  Sj  in  (3)  also  only  need  to  be  evaluated  at  the  integrator  step  times. 
Values  of  the  osculating  elements  a,  at  request  times  not  coinciding  with  the  integrator  step 
times  can  be  computed  from  (1)  using  interpolation  formulas. 

In  subsequent  chapters  we  will  outline  the  methods  of  derivation  and  give  explicit  formulas 
for  the  terms  A,,  C] ,  5/  corresponding  to  various  perturbing  forces. 
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2  Mathematical  Preliminaries 

2.1  Equinoctial  Elements 

The  generalized  method  of  averaging  can  be  applied  to  a  wide  variety  of  orbit  element 
sets.  The  equinoctial  elements  were  chosen  for  SST  because  the  variational  equations  for 
the  equinoctial  elements  are  nonsingular  for  all  orbits  for  which  the  generalized  method  of 
averaging  is  appropriate-namely,  all  elliptical  orbits. 

In  this  chapter  we  give  an  overview  of  the  equinoctial  elements,  which  are  osculating 
(even  though  they  do  not  have  hats).  They  are  discussed  in  more  detail  in  [Broucke  and 
Cefola,  1972],  [Cefola,  Long,  and  Holloway,  1974],  [Long,  McClain,  and  Cefola,  1975],  [Cefola 
and  Broucke,  1975],  [McClain,  1977  and  1978],  and  [Shaver,  1980]. 


2.1.1  Definition  of  the  Equinoctial  Elements 

There  are  six  elements  in  the  equinoctial  element  set: 


a  i 
a2 

03 

a4 

as 

«6 


a 

h 

k 

P 

Q 

\ 


semi  major  axis 

components  of  the  eccentricity  vector 

components  of  the  ascending  node  vector 
mean  longitude 


The  semimajor  axis  a  is  the  same  as  the  Keplerian  semimajor  axis.  The  eccentricity  vector 
has  a  magnitude  equal  to  the  eccentricity  and  it  points  from  the  central  body  to  perigee. 
Elements  h  and  k  are  the  g  and  f  components,  respectively,  of  the  eccentricity  vector  in 
the  equinoctial  reference  frame  defined  below.  The  ascending  node  vector  has  a  magnitude 
which  depends  on  the  inclination  and  points  from  the  central  body  to  the  ascending  node. 
Elements  p  and  q  are  the  g  and  f  components,  respectively,  of  the  ascending  node  vector  in 
the  equinoctial  reference  frame. 

There  are  actually  two  equinoctial  element  sets:  the  direct  set  and  the  retrograde  set.  As 
the  names  imply,  the  direct  set  is  more  appropriate  for  direct  satellites  and  the  retrograde  set 
is  more  appropriate  for  retrograde  satellites.  It  is  possible,  however,  to  use  direct  elements 
for  retrograde  satellites  and  vice  versa,  and  for  non-equatorial  satellites  this  presents  no 
problem.  For  equatorial  satellites  there  are  singularities  which  must  be  avoided  by  choosing 
the  appropriate  equinoctial  element  set.  For  direct  elements 

lim  Jp2  +  ?2  =  oo  (1) 

i— nr  v 

while  for  retrograde  elements 

Jim  yjp*  +  q2  =  oc  (2) 

For  each  equinoctial  element  set  there  are  three  associated  vectors  (f,  g,  w)  which  define 
the  equinoctial  reference  frame.  These  vectors  form  a  right-handed  orthonormal  triad  with 
the  following  properties: 
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Figure  2:  Retrograde  Equinoctial  Reference  Frame 


1.  f  and  g  lie  in  the  satellite  orbit  plan<\ 

2.  w  is  parallel  to  the  angular  mr'-  .entum  vector  of  the  satellite. 

3.  The  angle  between  f  and  the  ascending  node  is  equal  to  the  longitude  of  the  ascending 
node. 

This  leaves  two  choices  for  f  and  g,  one  associated  with  the  direct  element  set  and  one 
associated  with  the  retrograde  element  set.  The  two  sets  of  (f,  g,  w)  are  illustrated  in 
Figures  1  and  2.  In  the  Figures,  and  throughout  this  document,  (x,y,z)  denote  a  set  of 
Cartesian  coordinates  whose  origin  moves  with  the  center  of  mass  of  the  central  body  and 
whose  axes  are  nonrotating  with  respect  to  inertial  space. 


2.1.2  Conversion  from  Keplerian  Elements  to  Equinoctial  Elements 

If  a,  e,  i,  fi,  u>,  M  denote  the  conventional  Keplerian  element  set  then  the  equinoctial  ele¬ 
ments  are  given  by: 

a  —  a 

h  =  esin{w  +  /fi) 
k  =  ecos(u  +  Id) 

p  =  [tan  ^5)]  sin  ST  ^ 


q  =  [tan  (t)]*  cos  ft 
A  —  A/  -f-  w  ■+■  Id 

The  quantity  1  is  called  the  retrograde  factor  and  has  two  possible  values: 


for  the  direct  equinoctial  elements 
for  the  retrograde  equinoctial  elements 


There  are  two  auxiliary  longitudes  associated  with  the  equinoctial  element  set:  the  ec¬ 
centric  longitude  F  and  the  true  longitude  L.  They  are  related  to  the  Keplerian  eccentric 
anomaly  E  and  true  anomaly  /  by  the  equations: 

F  —  E  tu  -f-  Id  (3) 

L  =  f  +  u +  1(1  (4) 


These  auxiliary  longitudes  are  used  in  converting  from  equinoctial  elements  to  position  and 
velocity.  In  addition,  certain  perturbations  are  modeled  with  Fourier  series  expansions  in  F 
or  L. 


2.1.3  Conversion  from  Equinoctial  Elements  to  Keplerian  Elements 

In  order  to  convert  from  equinoctial  to  Keplerian  elements,  it  is  first  necessary  to  compute 
an  auxiliary  angle  which  is  defined  by: 


sin(  = 


cos  (  = 


fh2  +  k 2 
k 

< h 2  +  it2 


The  Keplerian  elements  are  then  given  by: 


I h 2  +  k 2 


=  t r  ^  ---)  +  2/ arctan  \Jp2  +  q2 


sin  fl  = 


cos  fi  = 


P 

s/P1  +  9* 
<7 

v/p2  +  9* 

=  <-/« 


M  =  A-C 

where  /  is  defined  by  (2. 1.2-2). 

The  Keplerian  eccentric  and  true  anomalies  are  given  by: 

£  =  F-C 
/  =  W 


2.1.4  Conversion  from  Equinoctial  Elements  to  Position  and  Velocity 

The  first  step  in  converting  from  equinoctial  elements  to  position  and  velocity  is  to  determine 
the  equinoctial  reference  frame  basis  vectors  (f,  g,  w).  Their  components  in  the  ( x ,  y,  z) 
system  are 

!  r  1  -  P2  +  <r  ' 

f  =  f  r>  ,--2  2p9 

1+p  [  ~2Ip 

i  r  2/p?  ] 


s  =  rTF77[(1%-,2)/ 

1  2p 

w  = - - - -  —  2q 

‘+P2  +  ?J 
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The  second  step  is  to  find  the  eccentric  and  true  longitudes  F  and  L ,  respectively.  To 
find  the  eccentric  longitude  F,  one  must  solve  the  equinoctial  form  of  Kepler’s  Equation  (see 
Section  6.1)): 

A  =  F  +  hcosF  —  ksinF  (2) 

Then  define  two  auxiliary  quantities,  the  Kepler  mean  motion  n  and  a  quantity  called  b: 


n 

b 


1 

1  +  yA  -  h2~-  k 2 


(3) 

(4) 


Here,  and  throughout  this  document,  n  is  the  gravitational  constant  GM  of  the  central  body. 
The  true  longitude  L  is  then  given  by: 


sinZ 
cos  L 


(1  —  k2b)  sin  F  +  hkb  cos  F  —  h 
1  —  h  sin  F  —  k  cos  F 
(1  —  h2b)  cos  F  +  hkb  sin  F  —  k 
1  —  h  sin  F  —  k  cos  F 


(5) 


The  third  step  is  to  compute  the  position  and  velocity  components  (X,  Y)  and  (X,  Y ) 
of  the  satellite  in  the  equinoctial  reference  frame.  The  radial  distance  of  the  satellite  is  given 
by: 

r  =  a(l  —  h  sin  F  —  fccosF)  =  ^ 


1  +  h  sin  L  +  k  cos  L 

The  position  components  are  then  given  by: 

X  =  a[(l  —  h2b)  cos  F  +  hkb  sin  F  —  k]  =  r  cos  L 
Y  =  a[(l  —  A:26)  sin  F  +  hkbcos  F  —  h]  =  r  sin  L 

The  velocity  components  are  then  given  by: 

v  na2r i r  n  l2l\  ■  fi  na(h  +  sin L) 

X  =  - [hkbcos  F  —  (1  —  h  b)  sin  F  = - .  ■  =4 

r  1  '  J  Vl~h2~  V 

T>  nfl2r/i  j.2i\  F  ...  •  ni  na(k  +  cos L) 

Y  =  - [(1  -  Jr6)cosF  -  hkbsmF  =  ~—===L 

r  u  ;  J  s/\-h2-  k2 

Here  dots  denote  differentiation  with  respect  to  the  time  t. 

The  final  step  is  now  to  compute  the  position  and  velocity  vectors: 


(6) 


(7) 


(8) 


r  =  Xf  +  Yg 
f  =  Xf  +  Kg 


(9) 
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2.1.5  Conversion  from  Position  and  Velocity  to  Equinoctial  Elements 

The  first  step  in  converting  from  position  r  and  velocity  r  is  to  compute  the  semimajor  axis 
a,  which  is  obtained  by  inverting  the  well-known  energy  integral  for  the  two  body  problem: 


a  = 


2_ 

Irl 


IfT 


(1) 


The  second  step  is  to  compute  the  basis  vectors  (f,  g,  w)  of  the  equinoctial  reference 
frame.  The  w  vector  is  obtained  by  normalizing  the  angular  momentum  vector: 


rxr 


Equinoctial  elements  p  and  q  are  then  given  by: 

wx 


1  +  Iwz 

Wy 

1  +  lWz 


(2) 


(3) 


Vectors  f  and  g  are  then  computed  using  elements  p  and  q  in  equations  (2. 1.4- la,  lb). 

The  third  step  is  to  compute  the  eccentricity-related  quantities.  The  eccentricity  vector 
e  is  given  by: 


e  = 


rx(rxr) 


(4) 


Equinoctial  elements  h  and  k  are  then  given  by: 


h  =  eg 
k  —  e  •  f 


(5) 


The  final  step  is  to  compute  the  mean  longitude  A.  First  compute  the  position  coordinates 
of  the  satellite  in  the  equinoctial  reference  frame: 

X  =  r  f 
Y  =  rg 

Then  compute  the  eccentric  longitude  F: 


sin  F 
cos  F 


L  (1  -  h2b)Y  -  hkbX 

+  ay/1  -  F  -  k2 

(1  -  k2b)X  -  hkbY 

+  ay/ 1  -  ft2  -  k* 


(7) 


where  b  is  defined  by  (2. 1.4-4).  The  mean  longitude  A  is  then  given  by  the  equinoctial  form 
(2. 1.4-2)  of  Kepler’s  equation. 
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2.1.6  Partial  Derivatives  of  Position  and  Velocity  with  Respect  to  the  Equinoc 
tial  Elements 

Let 


A  =  na2  =  y/jla 


1 


B  =  Vl  -  h2  -  k*  =  ^  -  1 
C  =  1  +  p2  +  ?2 


(1) 


The  partial  derivatives  of  the  position  vector  r  with  respect  to  the  equinoctial  elements  are 
then  given  by: 

r 


where 


dr 

da 

dr 

dh 

dr 

dk 

dr 

dp 

dr 

dq 

dr 

dX 

dX 


a 

dXe  dY 
■f+^ 


dh 

dX 


,  dY 
dkf+  dkg 
2[Iq(Yf-Xg)-Xv,} 


2I\p(Xg  -  Y f)  +  Yw] 
C 
r 
n 


(2) 


kX 


aYY 


~dh 

n{l  +  B) 

+  AB 

dY 

kY 

aXY 

dh 

n(l  +  B) 

AB 

dX 

hX 

—  _  j 

aYX 

dk 

n(l  4-  B)  X 

AB 

dY 

KY 

aXX 

dk 

n(l  +  B) 

AB 

—  a 


(3) 


The  partial  derivatives  of  the  velocity  vector  r  with  respect  to  the  equinoctial  elements 
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are  given  by: 


dX.  dY 
dhf+  dhg 
dXc  dY 
dk  f  +  dk  8 
2[lg(Yi  -  Xg)  -  Xw] 
C 

2I\p(Xg  -  Yf)  +  yw] 


where 


aY 2  A  /  akX  V2\ 

AB  +  r3  [l  +  B  ~  B  ) 
aXY  A  (  akY 
~  AB  +  r3  Vl  +  5+  B  ) 
aXY  A  (  ahX  XY\ 
AB  r3  \l  +  R  +  B  ) 
aX2  A  f  ahY  X2\ 
AB  r3Vl  +  5  BJ 


2.1.7  Partial  Derivatives  of  Equinoctial  Elements  with  Respect  to  Position  and 
Velocity 

Let  ^  and  represent  the  vectors  whose  components  in  the  (z,  y,  z)  system  are  the  partial 
derivatives  of  the  element  a,  with  respect  to  (z,  y,  z)  and  ( x,y,z ),  respectively; 

da j  da{  5a,  5a, 

5r  dxdydz 

5aj  _  5a,  5a,  5a, 

5r  dx  dy  dz 
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The  partial  derivatives  of  the  equinoctial  elements  with  respect  to  position  are  then  given 
by 


da 

2a2r 

dr 

r3 

dh 

abhBr 

l 

k(pX  ~  IqY)w 

B  dr 

dr 

3  + 

rJ 

AB 

A  dk 

dk 

dr 

= 

abkBr 

r3 

h{pX~IqY)w 

AB 

Bdr 

Adh 

dp 

CYw 

dr 

2  AB 

dq 

CX  w 

dr 

2AB 

d\ 

r  (PX 
A 

-  IqY)w  bB  f 

dr  .  dr 
hdh  +  kdk 

dr 

AB  A  l 

The  partial  derivatives  of  the  equinoctial  elements  with  respect  to  velocity  are  given  by: 


da  __  2r 
dr  A 

dh  (2 XY  -  XY)t  -  XXg  k(lqY  -  pX) w 

dir  p  +  AB 

dk  (2 XY  -  XY)g  -  YYf  h(IqY  -  pX )w 

dr  p  AB 

dp  CY  w 

di  ~  2 AB 

dq  __  ICX  w 

dir  ~  2AB 

flA  2r  tj-Af  (IqY-pX)  w 

dr  yl  1+  B  A 

2.1.8  Poisson  Brackets 


(3) 


The  Poisson  brackets  of  the  element  set  (<*!,•••,  ae)  are  defined  by  the  equations 

dat  daj  da,  da3 
(a„  a3)  — 

It  is  immediately  evident  that 


(a,,  ai )  =  0 

(an  aj)  =  ~{aji  ai) 


(1) 

(2) 
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The  fifteen  independent  Poisson  brackets  for  the  equinoctial  element  set  (a,  h,  k,  p,  q,  A) 
are  given  by 


(a,  h)  =  0 

<*•  A)  -  A(l  +  B) 

(a,  k)  =  0 

<*’  P)  =  2 AB 

(a,  p)  =  0 

>  hqC 
(  2AB 

(a,  q)  =  0 

,  kB 

,  2 

#  ,  c2 

(°»  A)  = 

na 

iAB1 

(A,  k)  = 

<P’  A)  =  2 AB 

(h  n\  kl>C 

{h'p)~  2  AB 

(h  rri  kqC 

{h'  q)  ~  2AB 

(<7,  ~  2 AB 

2.1.9  Direction  Cosines  (a,/?, 7) 

The  conservative  perturbations  are  more  conveniently  described  by  the  direction  cosines 
(a,  0, 7)  of  the  symmetry  axis  rather  than  the  equinoctial  elements  p  and  q.  For  central- 
body  gravitational  spherical  harmonics,  let  zB  be  the  unit  vector  from  the  center  of  mass 
to  the  geographic  north  pole  of  the  central-body.  For  third-body  point  mass  effects  and 
shadowless  solar  radiation  pressure,  let  z B  be  the  unit  vector  from  the  center  of  mass  to 
the  third-body.  The  direction  cosines  of  z B  with  respect  to  the  equinoctial  reference  frame 
are  then  given  by 

a  =  zBf 

0  =  zjs-g  (1) 

7  =  ZB  •  W 

The  quantities  (a,0, 7)  are  not  independent  but  related  by  the  equation 

a2  +  02  +  72  =  1  (2) 

Note  that  (a,  0, 7)  are  functions  of  p  and  g,  since  the  unit  vectors  (f,  g,  w)  are  functions 
of  p  and  q  through  equations  (2.1. 4-1).  Note  also  that  (a,/?, 7)  are  functions  of  f,  since 
z B  is  a  slowly  varying  function  of  time.  If  the  vector  zB  along  the  geographic  axis  of  the 
central-body  is  parallel  at  epoch  to  the  2-axis  of  Figures  1  and  2,  then  the  direction  cosines 
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r 


of  z b  are  at  epoch 


a  = 


0  = 


7  = 


2  Ip 


1  +  p2  +  g2 

2g 

1  +  P2  +  92 
(1-P2-?2)/ 


1  +  p2  +  g2 

The  partial  derivatives  of  (a,/?, 7)  with  respect  to  p  and  q  are: 


da 

2{lqP  +  7) 

dp 

C 

da 

2/p/? 

dq 

C 

W 

2  Iqa 

dp 

C 

00 

2I(pa  —  7) 

dq 

C 

d~i 

2a 

dp 

C 

d~f 

210 

dq 

C 

(3) 


(4) 


2.2  Variation-of- Parameters  (VOP)  Equations  of  Motion 

The  Cartesian  equations  of  motion  for  an  artificial  satellite  in  an  inertial  coordinate  system 
are  [Battin,  1987]: 

f  =  ~j^  +  0  +  ™  (1) 

Here  r  is  the  position  vector  from  the  center  of  mass  of  the  central  body  to  the  satellite, 
r  =  3pr  is  the  acceleration  vector,  p  =  GM  is  the  gravitational  constant  of  the  central 
body,  q  is  the  acceleration  due  to  nonconservative  perturbing  forces  (atmospheric  drag, 
solar  radiation  pressure),  and  Ti  is  a  potential-like  function  called  the  disturbing  function 
from  which  one  can  derive  the  acceleration  due  to  conservative  perturbing  forces  (central- 
body  spherical  harmonics,  third-body  point-mass).  If  m  and  II  are  the  mass  and  potential 
energy,  respectively,  of  the  satellite,  then  the  disturbing  function  ft  is  given  by: 

_  n 

71  = - 

m 

In  order  to  apply  the  generalized  method  of  averaging,  it  is  necessary  to  convert  the 
equations  of  motion  into  a  form  giving  the  rates  of  change  of  the  satellite  orbit  elements 
as  a  function  of  the  orbit  elements  themselves.  The  equations  of  motion  resulting  from 
this  conversion  are  called  the  Variation-of-Parameters  (VOP)  equations  of  motion.  The 
derivation  of  these  equations  is  discussed  in  some  detail  in  [Cefola,  Long,  and  Holloway, 
1974]  and  [McClain,  1977]. 
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In  this  section  we  let  (ai,---a6)  =  (a,  h,  k,  p,  q,  A)  denote  the  osculating  equinoctial 
elements  (even  though  they  do  not  have  hats).  Then  the  VOP  equations  of  motion  turn  out 
to  be 

3a,  JU  dU 

ai  =  n6t6  +  —q-2Jai,aj)j^-  (3) 

Here  f  =  ^  is  the  velocity  vector,  n  —  is  the  Kepler  mean  motion,  and  6,6  is  the 
Kronecker  delta.  The  partial  derivatives  3a, /dr  are  given  by  equations  (2. 1.7-3),  and  the 
Poisson  brackets  (a,,  a j)  are  given  by  equations  (2. 1.8-2, 3). 

The  VOP  equations  of  motion  (3)  include  three  contributions  to  the  orbit  element  rates 
of  change.  The  two-body  part  is: 

a,  =  nS, *  (4) 

The  Gaussian  or  nonconservative  part  is: 


3aj 

a‘  =  W  * 


The  Lagrangian  or  conservative  part  is: 


.  _  A  /  *** 

a,  _  !>•,  a))da. 


(5) 


(6) 


In  the  remainder  of  this  document  it  will  be  convenient  to  discuss  these  contributions  sepa¬ 
rately,  but  they  must  be  added  together  to  obtain  the  total  orbit  element  rates  of  change. 

The  Lagrangian  part  of  the  VOP  equations  of  motion  contains  the  partial  derivatives  of 
the  disturbing  function  71  with  respect  to  p  and  q.  The  perturbations  which  contribute  to  7 Z 
are  not  conveniently  described  in  terms  of  p  and  q,  however.  For  these  functions,  it  is  better 
to  write  71  as  a  function  of  (a,  h ,  k,  A)  and  the  direction  cosines  (o,  0,  7)  of  the  symmetry 
axis  of  the  perturbation.  The  partial  derivatives  of  the  disturbing  function  7?.  with  respect 
to  p  and  q  can  then  be  obtained  by  applying  the  Chain  Rule: 

&R 

dp 

&R  _  dU  da  3713/3 
dq  da  dq  +  d0  dq  +  d'j  dq 

The  partial  derivatives  of  (a,  0,  7)  with  respect  to  p  and  q  are  given  by  (2. 1.9-4).  To  simplify 
the  notation,  let  us  again  use  the  auxiliary  quantities  A,  B,  C  defined  by  (2. 1.6-1).  Also,  let 
us  define  the  cross-derivative  operator 


dTl  da  dTl  d0  dTl  dj 

"F  r.  /i  n  4" 


da  dp 
dIZ  da 


d0  dp 


37  dp 
371 37 


(7) 


^  dll  adU 

n,a0-ad0  0dQ 


(8) 


Note  that  7 l,ap  =  —7 l,pa.  Then  the  partial  derivatives  of  71  with  respect  to  p  and  q  turn  out 
to  be 

dfc  2  /  _  _  \ 

dp  ~  C  +Iqll'a0  ) 


dp 

m 

dq 


2/ 

C 


(tI,i3i  +v7l  1O0  ^ 


(9) 
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With  this  notation,  equations  (6)  for  the  Lagrangian  part  of  the  VOP  equations  of  motion 
become: 


h  = 
k  = 

P  = 
<7  = 
A  = 


2a  dll 
A  dX 
B&R 


+ 


A  dk  AB 
WBdfc  h 


ay  I QHiQ-j  ^ 


\B(fK  h  (  _  _  _  \ 

.Adh  +  AB  V>n'°'1  ) 


hB  m 
>1(1  +  B )  dX 
kB  &R] 
+  A(l  +  B)  d\. 


C 

2  AB 

C 

2  AB 
2a  dll 


P  -H 

q  (n,hk  —H 


+ 


B 


A  da  >1(1  +  B) 


H 

)-n 

/,  &r  ,dn\  i  >  ^  F  _  \ 

[h  dh  +  k  dk )  +  ab  \p7Z'^  /<?7^  ) 


&R 
■  dX 
&R 
dX 


ya-l 


(10) 


2.3  Equations  of  Averaging 

The  Generalized  Method  of  Averaging  may  be  used  to  divide  the  VOP  equations  of  motion 
(2.2-3)  into  a  short-periodic  part  which  can  be  integrated  analytically  and  a  slowly- varying 
part  which  can  be  integrated  numerically  using  time  steps  several  orders  of  magnitude  longer 
than  the  time  steps  appropriate  for  integrating  the  untransformed  equations  of  motion.  The 
Generalized  Method  of  Averaging  and  other  perturbation  techniques  are  discussed  in  [Nayfeh, 
1973].  Only  a  summary  of  the  application  of  this  procedure  to  (2.2-3)  will  be  given  here. 
More  details  can  be  found  in  [Cefola,  Long,  and  Holloway,  1974],  [McClain,  1977],  [McClain, 
Long,  and  Early,  1978]  and  [Green,  1979]. 

To  apply  the  Generalized  Method  of  Averaging  we  first  assume  that  the  osculating  orbit 
elements  a,  are  related  to  a  set  of  mean  elements  a,  by  a  near-identity  transformation: 


OO 

a,  =  a,  +  ^2^Vi(a,h,k,p,q,X,t) 
i= i 


(1) 


Here  again,  the  indexed  variables  (ai,-*-,a6)  refer  to  the  equinoctial  orbit  elements 
(a,h,k,p,q,X)  and  hats  distinguish  the  osculating  elements  from  the  mean  elements.  The 
quantity  represents  a  small  short-periodic  variation  of  order  j  in  element  i.  The  quantity 
e  is  called  the  “small  parameter”  and  plays  the  role  of  a  variational  parameter  in  deriving  the 
Equations  of  Averaging.  (Note  that  the  superscript  j  is  used  in  the  symbol  to  designate 
a  power  and  in  the  symbol  to  designate  an  index.) 
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The  short-periodic  variations  are  ac turned  to  contain  all  of  the  high-frequency  components 
in  the  osculating  elements  dj,  so  that  the  mean  elements  a,  vary  slowly  with  time.  This 
requirement  can  be  expressed  by  the  following  two  sets  of  inequalities: 


1  da 
n  dt 


1  \dai 


n 


dt 


1_ 

n 


d\ 

dt 


n 


a 

1  for  i  =  2, 3,4,5 

1 


(2) 


where  n  is  the  Kepler  mean  motion,  and 

|dfc+1a I 


A*+1 

A  *+1 


dtk+l 

dk+lat 


dtk+ 1 


<  a 


<C  1  for  i  =  2, 3, 4, 5, 6 


(3) 


where  k  is  the  order  and  A  is  the  step  size  of  the  numerical  integrator.  Inequalities  (2)  ensure 
that  second-order  effects  will  be  small,  while  inequalities  (3)  ensure  that  the  integrator  errors 
will  be  small. 

Using  the  variational  parameter  c,  we  can  write  the  osculating  Cartesian  equations  of 
motion  (2.2-1)  as 


d2r 
dt 2 


"j*F  +  e(q+v*) 


(4) 


As  c  increases  from  0  to  1  the  resulting  motion  varies  smoothly  from  two-body  motion  to 
the  actual  motion.  The  osculating  VOP  equations  of  motion  (2.2-3)  then  become 


ddi 

dt 


n(a)6ie  -f  c 


dat  «  .  .  .&R, 

aTq-I>’a'>a5- 


3=1 


which  can  be  written  in  the  form 


ddi 

~dt 


n(d)6i6  +  e^(a,  h ,  k,p,  q,  A,  t ) 


(5) 

(6) 


The  terms  eF;  give  the  osculating  element  rates  of  change  due  to  the  perturbing  forces  as 
functions  of  the  osculating  elements. 

We  assume  the  following  form  for  the  mean  VOP  equations  of  motion: 


dai 

dt 


n{a)Si6  +  £  t}A\ (a,  h,k,p,q,t) 

3=1 


(7) 


The  terms  eJA2  give  the  mean  element  rates  of  change  due  to  the  perturbing  forces  as 
functions  of  the  mean  elements.  For  most  perturbations,  the  A 2  are  independent  of  the  fast 
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variable  A  as  indicated  in  (7).  For  central-body  resonant  tesseral  harmonics,  the  A\  are  also 
slowly- varying  functions  of  j A  —  mO  (see  section  3.3). 

The  osculating  rate  functions  Fi  and  the  mean  rate  functions  A*  may  be  explicit  functions 
of  the  time  t  because  the  perturbing  forces  may  change  with  time  when  the  satellite  position 
and  velocity  are  held  constant,  e.  g.,  due  to  motion  of  the  Moon. 

We  now  expand  the  osculating  rate  functions  in  a  variational  power  series  about  the  mean 
elements: 


Here 


Fi(a,h,k,p,q,\,t)  =  Fi(a,h,k,p,q,\,t)  +  £V fi(a,h,k,p,q,  A,t) 

j=i 


(8) 

(9) 

(10) 


Similarly,  we  expand  the  osculating  Kepler  mean  motion  about  the  mean  semimajor  axis: 


Here 


n(a)  =  n(a )  +  £  e3N3(a) 

>=i 


Nl 


N2 

N 3 


3  *7i  /  \ 


_3r7?  15  (r/l)2' 

2  a  8  a2 


3»i?  155M 

2a  4  a2 


35  (’ll)3 

16  a3 


n(o) 


(11) 


(12) 

(13) 

(14) 


Having  all  the  necessary  expansions,  we  can  now  derive  the  Equations  of  Averaging. 
First,  differentiate  (1)  with  respect  to  t  and  use  (7)  to  obtain  one  expression  for  the  osculating 
element  rates.  Next,  expand  the  functions  on  the  right  side  of  (6)  using  (8)  through  (14) 
to  obtain  another  expression  for  the  osculating  rates.  Then  equate  the  two  expansions  and 
require  that  they  be  equal  for  all  values  of  e  between  0  and  1.  Since  the  powers  of  c  are 
linearly  independent,  the  coefficients  of  t3  must  be  equal.  Taking  j  =  1,2,3, •••  yields  the 
Equations  of  Averaging  of  order  1, 2, 3,  •  •  •  respectively: 

A]  +  -^j-n(a)  4-  =  Fi(a,h,k,p,q,\,t)  +  N'S#  (15) 
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(16) 

(17) 


In  the  Equations  of  Averaging  shown  above,  the  osculating  rate  functions  F,,  the  mean 
rate  functions  Af,  and  the  short-periodic  variations  rjj  contain  effects  due  to  many  pertur¬ 
bations.  In  order  to  obtain  practical  expressions  for  A\  and  ijj,  it  is  convenient  to  partition 
the  Equations  of  Averaging  to  separate  the  effects  of  different  perturbations. 

The  first  step  in  partitioning  the  Equations  of  Averaging  is  to  express  the  osculating  rate 
functions  Fj  in  terms  of  the  contributions  Fto  of  the  separate  perturbations: 

*Fi  =  Yl/°Fia  (18) 

a 

The  sum  is  taken  over  all  perturbations  of  interest.  The  parameters  va  are  variational 
parameters,  one  for  each  perturbation.  Each  va  can  vary  from  0,  at  which  the  perturbation 
is  turned  off,  to  1,  at  which  the  perturbation  has  its  actual  strength.  The  partitioned 
Equations  of  Averaging  are  required  to  be  valid  for  all  values  of  the  i/Q. 

Substituting  equation  (18)  into  the  first-order  Equations  of  Averaging  (15)  leads  to  the 
following  expressions  for  A]  and  tj]  : 


eAj 

—  )  '  ^oAia 

o 

(19) 

= 

(20) 

a 


Substituting  equations  (18)— (20)  into  the  second-order  Equations  of  Averaging  (16)  leads  to 
the  following  expressions  for  Af  and  r/f: 

<?A2t  =  (21) 

0(3 

=  J2  WMioP  (22) 

o/3 

Substituting  equations  (18)— (22)  into  the  third-order  Equations  of  Averaging  (17)  leads  to 
the  following  expressions  for  Af  and  rjf: 

=  Y  v°v0v‘JAialh  (23) 

a(h 

=  Y,  Wpi'-iVioP-,  (24) 

0/37 

Similar  expressions  exist  at  higher  orders.  (Remember  that  the  first  index  in  Fia ,  Aia,  t?io, 
etc.  refers  to  the  orbit  element.) 
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If  we  substitute  (l8)-(20)  into  (15)  and  then  equate  the  coefficients  of  each  va,  we  obtain 
the  partitioned  form  of  the  first-order  Equations  of  Averaging.  Similar  procedures  with  (16), 
(17),  •  •  •  lead  to  higher-order  equations.  The  partitioned  Equations  of  Averaging  of  arbitrary 
order  cam  be  written  in  the  concise  form 


A,  +  §£»(«)  +  ^  =  Gi  - 

Explicit  expressions  for  the  Gi  functions  up  to  order  3  are: 

Ci^o  =  Fia(a,  h,  k,p,  q,  A,  t) 


(25) 


(26) 


Gic,0  =  53 


dJk 

da } 


Vjff  + 


15  l]la^}l0 
T  a2  " 


(«)*  6 


6 


E 


dyta 

da, 


A}0 


(27) 


Giap-y  — 


6 


e 


dF„ 

das 


i  6  6 

Vj0i  +  2  EE 

£  j=i  l 


d2F,„ 


+ 


f  15  ij iai]\0^ 
U  a2 


+ 


35 

16 


VlctV\0^l\'y 


(28) 


Comparing  the  partitioned  Equations  of  Averaging  (25)-(28)  with  the  full  Equations  of 
Averaging  (15)-(17),  we  see  that  the  first-order  equations  are  identical.  The  second-order  and 
higher-order  partitioned  equations  include  auto-coupling  equations  (e.g.,  a  =  ft  —  •  *  *  —  1) 
which  are  identical  to  the  full  equations,  but  in  addition  include  cross-coupling  equations 


(e.g.,  a  -  1,/?  =  2). 

The  partitioned  Equations  of  Averaging  (25)-(28)  give  the  fundamental  relations  which 
can  be  used  to  derive  expressions  for  the  mean  element  rates  A,q,  Ata0,  Axa0-,  and  the  short- 
periodic  variations  T^a,  t Then  the  total  mean  element  rates  Aj,  A,2,  Af  and  short- 
periodic  variations  Tjf ,  172,  rjf  can  be  obtained  from  the  decomposition  relations  (19)-(24). 


2.4  Averaged  Equations  of  Motion 

The  Equations  of  Averaging  (2.3-25)  can  be  solved  for  the  mean  element  rates 
Aj<*,  Aiapx  Aiafa  ■  •  •  by  applying  an  averaging  operator  <•••>,  to  be  defined  in  this  sec¬ 
tion,  to  both  sides  of  each  equation.  The  resulting  expressions  for  the  A, a,  Aia/J,  A,a/J-y 
can  then  be  added  together  as  shown  in  equations  (2.3-19,  21,  23)  and  substituted  into 
equations  (2.3-7)  with  e  =  va  =  1  to  form  the  mean,  or  averaged,  equations  of  motion: 

—■  =  n(a)6ft  +  53  +  53  +  53  Aar0Tr  +  •  ■  •  ( 1 ) 

O  a  0  o,0,-r 

These  equations  can  then  be  integrated  with  a  numerical  integrator  to  obtain  values  for  the 
mean  elements  ai  at  a  given  time. 
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The  averaging  operator  is  required  to  be  linear;  that  is,  if  p  and  a  are  any  two  real 
numbers  and  /  and  g  are  any  two  real  piecewise  continuous  functions  of  the  mean  elements, 
then: 

<  pf  +  ag  >=  p  <  f  >  +cr  <  g  >  (2) 

The  averaging  operator  is  also  required  to  be  idempotent;  that  is,  if  /  is  any  real  piecewise 
continuous  function  of  the  mean  elements,  then: 

«/>>  =  </>  (3) 

To  make  the  averaging  transformation  useful,  we  also  require  the  averaging  operator  to  have 
the  property  that  solving  the  Equations  of  Averaging  with  it  yields  slowly- varying  mean 
element  rates  and  small  short- periodic  variations. 

In  order  to  be  able  to  divide  the  Equations  of  Averaging  into  separate  equations  for  the 
mean  element  rates  and  the  short- periodic  variations,  we  impose  the  following  conditions: 

<tj,>=0  (4) 

It  is  not  immediately  obvious  from  the  Equations  of  Averaging  (2.3-25)  that  the  short- 
periodic  variations  can  be  required  to  average  to  zero  (equations  (4)).  Let  us  first  observe 
from  (2.3-26,  27,  28)  that,  at  any  order,  G,  are  predetermined  functions  of  the  osculating 
rate  functions  Fta  and  the  solutions  of  the  lower-order  Equations  of  Averaging.  Therefore  G, 
are  fixed,  while  we  can  vary  A,  and  77,  in  any  manner  which  satisfies  (2.3-25).  Let  us  assume 
that  the  short- periodic  variations  77,  do  not  average  to  zero,  and  write 

<  Vi  >=  (5) 


Let  us  then  define 


dki  dk,  3  n 

-  /t'  +  "aX+ aT  + 

=  Vi  ~ 


(6) 

(7) 


Solving  (6)-(7)  for  A,  and  77,  and  substituting  the  resulting  expressions  into  (2.3-25)  yields 

A'^ndl1'  4.  dr}*  r  3r?L/;  m 

A'  +  nd\  +  ~dt  G,~  27  6,6  (8) 

We  see  that  A\  and  77-  are  solutions  to  the  Equations  of  Averaging,  and  we  can  thus  choose 
A\  and  77,-  to  be  the  preferred  solutions.  Averaging  (7)  and  applying  (2,  3,  5)  yields 

<  V[  >=  0  (9) 


We  can  therefore  require  the  short-periodic  variations  to  average  to  zero  (equations  (4)). 

If  the  osculating  rate  functions  Fia  for  a  given  perturbation  are  small,  27r-periodic  in  the 
satellite  mean  longitude  A,  and  slowly- varying  in  time  when  the  satellite  orbit  elements  are 
held  fixed,  then  the  single-averaging  operator  has  the  required  properties: 

<  /  >=  ^  J  *f{a,h,k,p,q,\,t)d\  (10) 

Most  of  the  perturbations  commonly  acting  on  a  satellite  can  be  single-averaged: 
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1.  Central-body  gravitational  zonal  harmonics. 

2.  Third-body  gravitational  point  mass  effects. 

3.  Atmospheric  drag. 

4.  Solar  radiation  pressure. 

Some  perturbations  are  quickly-varying  when  expressed  as  a  function  of  time  but  slowly- 
varying  when  expressed  as  a  function  of  both  time  and  a  perturbing-body  phase  angle  0 
which  varies  linearly  with  time.  If  the  osculating  rate  functions  Fia  for  such  a  perturbation 
are  small,  2?r-periodic  in  both  A  and  9 ,  and  slowly-varying  in  time  when  6  and  the  satellite 
orbit  elements  are  held  fixed,  then  the  double- averaging  operator  has  the  required  properties: 


</>  =  f(aih,k,p,q,  A,  0,t)d\d0 


(11) 


4?r2  j —it  j—it 

+  2^2  51  [«*(;' A  ~ m  *>/:/:  /(a,  h,  k,p ,  q ,  A \9',t)  cos(j’A'  —  m9')d\'d9' 

(j,m)£B 


+  sin(jA  —  m9)  J  J  J(a,h,k>p,q,\' ,9' —  mO')d\'dO' 


(Equation  (11)  may  be  written  in  alternate  forms,  one  of  which  is  used  in  (3.3-2).)  Here  B 
is  the  set  of  all  ordered  pairs  (j,  m)  with  the  following  properties: 


m  >  1  (12) 

\j\  -  m0\  <  —  (13) 

T 

where  r  is  the  minimum  period  desired  for  perturbations  included  in  the  averaged  equations 
of  motion.  The  same  minimum  period  should  be  used  for  all  double-averaged  perturbations. 
Inequality  (13)  is  called  the  resonance  condition  and  denotes  the  satellite  frequencies  j  and 
perturbing-body  frequencies  m  which  are  in  resonance  with  each  other.  It  often  happens 
that  the  satellite  has  no  resonances  with  the  perturbing  body,  in  which  case  set  B  is  empty. 
The  minimum  period  r  should  obey  the  following  inequalities: 


r  >  3ta  (14) 

r  >  3 t8  (15) 


where  t\  and  Tg  are  the  periods  of  the  satellite  and  the  perturbing-body  phase  angle  0 
respectively.  In  addition,  r  should  be  at  least  8  times  as  long  as  the  step  size  used  in  the 
numerical  integrator,  and  much  longer  if  accurate  integration  of  perturbations  with  this 
period  is  desired.  It  is  dangerous  to  make  r  too  long,  however.  Components  of  the  double- 
averaged  perturbation  which  have  periods  equal  to  r  or  shorter  will  be  excluded  from  the 
averaged  equations  of  motion  and  treated  as  short-periodic  variations.  If  r  is  too  long  and 
deep  resonance  exists,  some  of  these  short-periodic  variations  may  be  large  enough  to  cause 
large  second-order  coupling  effects,  making  the  averaging  expansions  (2.3-1)  and  (2.3-7) 
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diverge.  For  this  reason,  r  should  be  less  than  100  times  the  integrator  step  size.  Making  r 
small  enough  to  ensure  convergence  of  the  averaging  expansions  takes  priority  over  inequality 
(15),  which  must  be  dropped  if  the  perturbing-body  phase  angle  0  varies  too  slowly.  If  0  is 
the  rotation  angle  of  the  Earth,  this  will  not  be  necessary,  since  the  Earth  rotates  quickly. 
If  0  is  the  rotation  angle  of  Mercury  or  Venus,  it  may  be  necessary  to  drop  condition  (15), 
depending  on  how  large  the  m-daily  (j  =  0)  short-periodic  variations  due  to  the  gravitational 
tesseral  harmonics  are.  The  Moon  is  a  borderline  case. 

The  following  perturbations  can  be  double-averaged: 

1.  Central-body  gravitational  sectoral  and  tesseral  harmonics.  For  this  perturbation,  0  is 
the  rotation  angle  of  the  central  body.  (For  a  more  precise  definition  of  0,  see  Section 
2.6.)  The  double-averaged  central-body  gravitational  spherical  harmonic  model  in  SST 
is  fully-developed  and  will  be  discussed  in  Section  3.3. 

2.  Third-body  point-mass  effects,  if  the  third  body  orbits  the  central  body.  For  this 
perturbation,  0  is  the  equinoctial  mean  longitude  of  the  third  body.  The  current  double- 
averaged  third-body  model  assumes  that  the  third  body  orbits  the  central  body  in  a 
slowly- varying  Keplerian  ellipse.  Methods  for  predicting  the  effects  of  short-periodic 
variations  in  the  third-body  orbit  on  the  satellite  orbit  have  yet  to  be  developed.  For 
Earth  satellites,  the  short-periodic  variations  in  the  orbit  of  the  Moon  are  substantial. 
If  they  are  included  in  the  Lunar  ephemeris  used  with  the  double- averaged  third-body 
model,  the  integrator  step  size  will  be  driven  down  to  values  appropriate  for  the  single- 
averaged  third-body  model,  thus  destroying  the  usefulness  of  the  double- averaging 
expansion.  The  step  size  reduction  is  avoided  by  using  a  smoothed  ephemeris  for  the 
Moon,  but  this  creates  an  error  in  the  computed  satellite  mean  element  rates  due  to 
the  Lunar  perturbation,  and  the  size  of  this  error  is  not  known.  Because  of  these 
limitations,  double-averaged  third-body  perturbation  models  will  not  be  discussed  in 
detail  in  this  document.  For  a  complete  description  of  the  current  model,  see  [Collins, 
1981], 

There  are  some  perturbations  for  which  no  averaging  operator  with  the  required  proper¬ 
ties  can  be  found.  These  perturbations  are  called  non-averageable  and  include: 

1.  Atmospheric  drag,  with  an  asymmetric  spacecraft  and  fast,  non-periodic  attitude  vari¬ 
ation. 

2.  Solar  radiation  pressure,  with  an  asymmetric  spacecraft  and  fast,  non-periodic  attitude 
variation. 

3.  Continuous  thrust,  with  fast,  non-periodic  changes  in  direction. 

4.  Impulsive  thrust. 

These  perturbations  are  typical  of  directed  flight,  which  in  general  is  not  required  to  be  either 
slowly-varying  or  27r-periodic  in  A  or  any  other  phase  angle.  Semianalytic  Satellite  Theory 
cannot  predict  the  effects  of  these  perturbations,  with  the  exception  of  impulsive  thrust  (see 
below),  and  should  not  be  used  unless  they  are  small  enough  to  ignore. 
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Some  types  of  directed  flight  are  averageable,  and  these  include  many  scenarios  of  prac¬ 
tical  interest: 

1.  A  drag- perturbed  satellite  whose  solar  panels  always  point  directly  at  the  Sun.  The 
attitude  of  the  satellite  relative  to  the  atmosphere  will  be  27r-periodic  in  A  and  will 
vary  slowly  in  time  as  the  Earth  revolves  about  the  Sun. 

2.  A  spacecraft  with  a  solar  sail  which  is  feathered  when  approaching  the  Sun  and  per¬ 
pendicular  to  the  Sun  line  when  receding.  The  attitude  of  the  sail  will  be  27r-periodic 
in  A  and  will  vary  slowly  in  time  as  the  Earth  revolves  about  the  Sun. 

3.  A  spacecraft  with  a  constant-thrust  ion  engine  whose  thrust  is  always  parallel  to  the 
orbit.  The  direction  of  the  resulting  acceleration  will  be  27r-periodic  in  A  and  the 
magnitude  will  vary  slowly  in  time  as  reaction  mass  is  depleted  and  the  spacecraft  gets 
lighter. 

Of  course,  the  perturbations  remain  averageable  only  as  long  as  the  orbit  remains  elliptical. 
Parabolic  and  hyperbolic  orbits  are  beyond  the  scope  of  this  document. 

Impulsive  thrust  is  not  averageable,  but  its  effects  can  be  predicted  using  the  following 
procedure: 

1.  Integrate  the  averaged  equations  of  motion  (1)  up  to  the  time  of  the  impulsive  maneu¬ 
ver. 

2.  Compute  the  short-periodic  variations  as  functions  of  the  mean  elements  using  equa¬ 
tions  (2.3-20,  22,  24)  and  the  equations  in  Section  2.5. 

3.  Add  the  short-periodic  variations  to  the  mean  elements  to  get  the  osculating  elements 
(equations  (2.3-1)). 

4.  Convert  the  osculating  equinoctial  orbit  elements  to  position  and  velocity  using  the 
equations  in  Section  2.1.4. 

5.  Add  the  velocity  change  Av  caused  by  the  impulsive  maneuver  to  the  satellite  velocity. 

6.  Convert  the  satellite  position  and  velocity  to  osculating  equinoctial  orbit  elements  using 
the  equations  in  Section  2.1.5. 

7.  Invert  equations  (2.3-1)  to  convert  the  osculating  elements  to  mean  elements  (see  Sec¬ 
tion  6). 

This  procedure  is  always  valid,  but  if  impulsive  thrusts  occur  more  often  than  once  per  orbit, 
it  may  be  too  expensive  to  make  the  use  of  averaging  worthwhile.  For  additional  methods 
for  modeling  impulsive  maneuvers,  as  well  as  continuous  thrust,  see  McClain  [1982]. 

Some  perturbations  are  usually  averageable,  but  cannot  be  averaged  in  certain  circum¬ 
stances  because  of  large  second-order  effects.  These  include: 
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1.  Third-body  point-mass  effects,  for  a  satellite  whose  orbit  comes  too  close  to  the  orbit 
of  the  third  body.  This  perturbation  will  be  non-averageable  even  if  resonance  locking 
ensures  that  the  satellite  will  always  remain  far  from  the  third  body. 


2.  Third-body  point-mass  effects,  for  a  satellite  whose  orbit  comes  too  close  to  the  bound¬ 
ary  of  the  central  body’s  gravitational  Sphere  of  Influence. 


3.  Atmospheric  drag  effects,  during  the  terminal  stage  of  reentry. 

Semianalytical  Satellite  Theory  should  not  be  used  under  these  circumstances. 

It  is  worth  considering  at  this  point  whether  the  averaging  operators  (10)  and  (11)  ac¬ 
tually  have  the  required  properties  (2)-(3).  It  is  immediately  clear  from  equations  (10)  and 
(11)  that  both  operators  are  linear  (equation  (2)).  It  is  also  clear  from  equation  (10)  that  the 
single-averaging  operator  is  idempotent  (equation  (3)).  To  show  that  the  double-averaging 
operator  is  idempotent,  we  double-average  both  sides  of  equation  (11)  to  obtain 


dX  dO 


«  f  »  =  T— 2  /  f  <  /  > 

4 7T2  J-i,  J-n 

+  ~r  Y  [cos(j\  —  m0)f  f  <  /  >  cos(jA'  —  m9')dX'dQ' 
2n  (j.m)66  Jr~\ 

-fsin(jA  —  mO)  J  J  </  >  sin(jA'  —  mQ')d\'dO'] 
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T.  [cos(jA  —  m0 )  f  I  cos2(jA' 
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-  mQ')dX'dO'  J*  j '  fcos(j\"  -  m6")d\"d6" 

-  m0')dX'd0'  J*  J W  /sinO'A"  -  m0")d\"d0"] 


(ifl) 

Note  that,  if  /  is  independent  of  0,  the  double- averaging  operator  (11)  reduces  to  the 
single-averaging  operator  (10).  Thus  we  can  without  inconsistencies  apply  the  single-averaging 
operator  (10)  to  perturbations  that  do  not  depend  upon  6,  and  the  double- averaging  operator 
(11)  to  perturbations  that  do  depend  upon  9. 

If  the  perturbations  can  be  averaged  with  the  operators  (10)  or  (11),  then  the  Equations 
of  Averaging  (2.3-25)  can  be  solved  for  the  mean  element  rates  by  averaging  both  sides  of 
each  equation.  Note  that  the  averages  are  always  taken  with  ( a,h,k,p,q,t )  held  constant. 
Using  (4),  we  thereby  obtain  for  the  mean  element  rates  A;  at  any  order 


Ai  =<  G,  >  (17) 

Substituting  (2.3-26,  27,  28)  into  (17),  we  can  obtain  explicit  expressions  for  the  mean 
element  rates  at  each  order. 

The  explicit  equations  for  the  first-order  mean  element  rates  have  the  same  form  for  both 
single-averaged  and  double- averaged  perturbations: 


Aia  =<  Fia(a,h,k,p,q,\,t )  > 


(18) 
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The  averaging  operator  used  is  single  in  the  first  case  and  double  in  the  second. 

The  explicit  equations  for  the  second-order  mean  element  rates  involve  coupling  terms 
between  two  perturbations,  which  may  be  two  different  perturbations  (cross- coupling)  or 
copies  of  the  same  perturbation  (auto-coupling).  If  either  perturbation  is  single- averaged, 
the  equations  take  the  form 

Axa0  =  +  "g"^2  (fl aVw)  (19) 

If  both  perturbations  are  single- averaged,  a  single- averaging  operator  is  used.  If  either 
perturbation  is  double- averaged,  a  double-averaging  operator  is  used.  If  both  perturbations 
are  double- averaged,  the  equations  take  the  form 

A-  =  £{%?»}  +  T?  6*  - 1  (%*»)  <*» 

If  both  perturbations  use  the  same  perturbing-body  phase  angle  6,  then  a  double-averaging 
operator  is  used.  However,  if  the  perturbations  use  different  perturbing-body  phase  angles 
0  and  6\  then  a  triple- averaging  operator  is  used. 

Derivation  of  the  explicit  equations  for  the  third-order  mean  element  rates  will  be  left  as 
an  exercise  for  the  reader. 


2.5  Short-Periodic  Variations 

Once  the  mean  element  rates  A,a,  Alap,  Aiap^ ,  •  •  •  are  known,  the  Equations  of  Averaging 
(2.3-25)  can  be  solved  for  the  short-periodic  variations  */,■„,  t/IO(9,  TjtQp^,  •  •  •  by  using  Fourier 
series  expansions  for  the  functions  Gto,GiO0,Glo/fj^,*  •  -.  The  resulting  expressions  for  the 
Vian  Viaffi  Via0~ri ' ' ‘  can  then  t>e  added  together  as  shown  in  equations  (2.3-20,  22,  24)  with 
c  =  va  =  1  and  substituted  into  equations  (2.3-1)  with  e  =  1  to  give  the  osculating  elements: 

hi  =  a,  +  ^2  Via  +  Viap  +  JZ  Via0y  H -  (1) 

°  O0  Q0-y 

Of  course,  using  Fourier  series  expansions  for  the  functions  G,  assumes  that  the  osculating 
rate  functions  Fia  are  27r-periodic  in  the  phase  angles  of  the  expansions.  Most  perturbations 
can  be  expressed  with  more  than  one  kind  of  Fourier  series  expansion.  In  this  section  we  shall 
give  formulas  for  the  short-periodic  variation5;  corresponding  to  several  possible  expansions. 

Some  of  the  general  expressions  presented  in  this  section  are  new,  so  enough  steps  in 
the  derivations  will  be  presented  to  enable  the  reader  to  rederive  them.  Our  derivations 
are  based  on  the  work  of  [Cefola  and  McClain,  1978],  [Green,  1979],  [McClain  and  Slutsky, 
1980],  [Slutsky,  1980],  [Slutsky  and  McClain,  1981],  and  [McClain,  1982]. 


2.5.1  General  rji  Expansions  in  A 

Recall  that  the  Equations  of  Averaging  of  arbitrary  order  can  be  written  as 


A  _l  drii 

‘4i  +  "aX 


drji  3  n 

~di  =  Gi~  2Z6xm 


(i) 
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where  G;  are  predetermined  functions  (equations  (2.3-26,  27,  28))  and 

At  =<  G,  >  (2) 

At  order  m,  the  Equations  of  Averaging  contain  m  perturbations,  which  may  all  be  different 
or  may  include  multiple  copies  of  the  same  perturbation. 

If  all  of  the  perturbations  in  the  Equations  of  Averaging  (1)  are  single- averaged,  then  the 
equations  can  be  solved  for  the  short-periodic  variations  t]i  by  integrating  over  the  satellite 
mean  longitude  A.  It  is  convenient  to  first  define  the  short-periodic  kernels 

(3) 

The  constant  of  integration  is  specified  by  requiring 

<  ti  >=  o  (4) 


In  the  absence  of  explicit  time-dependence,  the  full  short-periodic  variations  rji  are  then  given 
by: 

/  £,  dx  (5) 

The  conditions  (2.4-4)  and  (4)  require 


(/<■")- 


0 


(6) 


In  the  presence  of  explicit  time-dependence,  the  full  short-periodic  variations  t /,  are  given 
by: 

K  (-l)fc  r&ti 


*  - 

-&4l('dX  +  gk  +  1)(^LlFdXW. 


(7) 


Here  K  is  the  order  of  the  highest  order  partial  derivatives  desired  in  the  expansions.  The 
conditions  (2.4-4)  and  (4)  require 


i^J  =  0  for  all  k  > 


1 


(8) 


Here  the  following  somewhat  unusual  notation  is  used  for  the  operator  which  is  the  kth 
indefinite  integral  (the  inverse  of  the  kth  derivative  operator): 


Jkf{X)dXk  =  J-J  f(X)dX  ■■■dX, 


(9) 


Alternately,  knowing  the  short-periodic  variations  T]f  in  the  absence  of  explicit  time- 
dependence,  we  can  compute  recursively  the  short-periodic  variations  r/f  including  the  kih 
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order  time  derivatives: 

1°  =  (i 

’•*  - + ^  L  ^ + L  i£dxM  <10) 

Explicitly,  we  suppose  that  the  functions  Gi  can  be  written  as  a  Fourier  series  in  the 
mean  longitude  A: 

Gi(a,h,k,p,q,\,t)  =  C?(a,h,  k,p,q,t)  (11) 

+  £  [C,-(a,h,  fc,p,  g,f)cos;A  +  S-(a,h,k,p,  q,t)sin  j'A] 
j=i 

where 

c'°  -  s/>" 

C\  =  -  T  G,  cos  jXdX  (12) 

j r  J— » 

5/  =  —  f  GisinjXdX 

7T  —  7T 

Using  (2.4-10)  and  (2),  we  can  obtain  the  mean  element  rates  A{-. 

Ai  =  C?  (13) 

We  can  then  obtain  the  short-periodic  variations  t/,  by  performing  the  integrals  (7,  8).  The 
final  result  of  these  calculations  is: 


r)i  =  £[C/  cos; A  -f  5-  sin; A] 
i- 1 


where 


JL\S, 

i  3 1U  1 

2a;  V(;n)2 

ac; 

a 

jn  ' 

1 

325/  3  3  <5,6  32C{ ' 

1 

(jn)3 

3<2  2a;  3f2 

(;'n)4 

5/  =  4-  C?  +  2i^5il  +  - 


2a  j 


1  [35/  32^6^1 


(;n)2  [  dt  2  a  j  dt 


1  32C,-  33^3^1  _  _J_  [335/  _  3 4 6# dPC{ 

(jn)3  3<2  2a  ;  3t2  (;n)4  3<3  2a  ;'  3f3 


2.5.2  General  »?,  Expansions  in  F 

In  this  section  we  suppose  that  the  functions  Gi  are  expanded  as  a  finite  Fourier  series  in 
the  eccentric  longitude  F : 

Gi{a,h,k,p,q,F,t)  =  C?(a,h,k,p,q,t)  (1) 
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where 


J 

+  YL  [c/(a> h'  P»  0  cos  jF  +  S;f  (a,  h,  fc, p,  $,  t )  sin  jF] 


<?  -  hL°<iF 

c\  =  -  r  Gi  cos  jFdF  (2) 

ir  j-it 

5/  =  -  I*  Gi  sin  jFdF 

7T  J 

Again  we  may  use  (2.4-10)  and  (2. 5. 1-2,  3,  4,  7,  8)  to  obtain  the  short  periodic  variations 
rjt.  To  do  this,  it  will  be  necessary  to  convert  the  integrals  over  A  into  Fourier  series  expansions 
over  F.  We  begin  by  supposing  /(A)  has  a  finite  Fourier  series  expansion  in  F  with  known 
coefficients  and  it  averages  to  zero: 

j 

/(A)  =  C°  +  cos  jF  +  S>  sin  jF)  (3) 

j= i 


<  /(A)  >=  0 


Using 


/  f(X)dX  =  J  m^dF 

the  following  consequence  of  the  equinoctial  form  (2. 1.4-2)  of  Kepler’s  Equation 

d\  r  .  _ 

— —  =  —  =  1  —  n  sin  F  —  k  cos  F 
oF  a 

and  the  following  well-known  trigonometric  identities 

_ ■  p  _  ,  r  «»(;'  -  k)F  +  cos(;  +  k)F 


cos  jiFcos  kF  = 
cos  jF  sin  kF  = 


sin(j  -f  k)F  —  sin(j  —  k)F 


„•  c> _ l  ip  _  sin(j  +  k)F  +  sin(;  -  k)F  ] 

sin  ?r  cos  let  —  - 

2 

.  .  „  .  ,  „  cos  (j  —  k)F  —  cos(j  +  k)F 

smiFsmkF  =  — - - 

2 

we  can  convert  the  right  side  of  (5)  into  a  Fourier  series  expansion  in  F.  Note  in  particular 
that  the  condition  (4)  implies  that  the  constant  term  C°  in  (3)  must  be  related  to  the  Fourier 
coefficients  C1  and  S1  by 

C°  =  k-C'  +  jS1  (8) 

Higher-order  integrals  can  be  computed  using  recursion  formulas  obtained  from  the  equa- 
tion 

/  /(A)<iA"tl  =  (9) 

«'m+l  J  Jm 
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We  summarize  the  conversion  for  the  general  multiple  integral  with  the  following  notation: 


/  f(\)d\m  =  Ul(C<,S<)  +  Jj2[UL(CC,S<)cosjF+VMC<,S<)smjF]  (10) 

Jm  >=i 

<  f  >=  0  (11) 

Jm 

Here  the  arguments  (C<‘,5I‘)  denote  that  the  functions  and  depend  upon  the  coeffi¬ 
cients  CJ  and  SJ  appearing  in  (3)  through  the  relations: 
for  1  <  j  <  J: 

vi=c> 

Vi  =  5J 

for  m  >  0  : 

K  =  + 

v'. ,  =  —u'-{i-—\v'--v2+-v2 

^m+i  2  ^ m  l  1  2  J  r m  2  m  '  2  v m 

<»  =  (i  -  j)  ul  -  yv~  ~ 

(12) 

for  2  <  j  <  J  +  m  +  1 : 

}  (-vi  +  fvi-'  +  jvs-  -  +  jvr1  j 

-  -CF"1  +  -1/J-1  -  -//J+1  _  -Vrj+I^ 
m+i  j  y  m  2  '  2  2 u  m  2  * m  J 

for  j  >  J  +  m  +  1: 


uL  =  o 

*2  =  0 


With  this  notation,  we  can  write  the  mean  element  rates  A,  and  the  short-periodic 
variations  j/<  corresponding  to  (1): 


A, 

»7. 


J+K+2 

=  £?+£)[<??  cos 
j- 1 


jF  +  S'/  sin;F] 


(13) 

(14) 
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where 


C\ 

si 


c° =  5C‘  +  5s- 


l-vi(^)+it^ww0,  2$) 

-&*  (<$•*) + |<* + fis&  ^)] 

nVl  ^ ’5* ) +  g -^rVk+1  V  aF’  IF j 

f-vJ  (cc  sc)  +  T(k  +  nizll^v"  (*£  *£\ 


_1, 

2a6*6 


(15) 


2.5.3  General  77,  Expansions  in  Z, 

In  this  section  we  suppose  that  the  functions  G,  are  expanded  as  a  finite  modified  Fourier 
series  in  the  true  longitude  L: 


Gt(a,h,k,p,  q,L,t ) 


M 

C?{a,h,k,p,q,t)  +  Y,  T>™{a,h,k,p,q,t){L-  \)m 

m=l 

J 

+  LP<«,M  ,p,q,t)cosjL  +  S-(a,h,k,p,  q,t)sin  jZ,] 
i= i 


(1) 


Here  the  quantities  (L  —  A)m  are  written  separately,  rather  than  by  replacing  them  with  their 
Fourier  series  expansions  recorded  below. 

The  equation  of  the  center  may  be  calculated  from  the  Fourier  series  expansion 


Here 


L  —  A  =  ^2  -(a j  cos  jL  —  p3  sin  jL) 
j=i  J 


Pi 


i  r 

—  (cos  jL)  =  —  cos  jLdX 
2ir  J —it 

1  r* 

(Tj  =  (sin  jL)  =  —  /  sin  jLdX 

27T  */— tt 

The  auxiliary  quantities  pj  and  (Tj  can  be  evaluated  using  the  equations  * 


‘Let  Z  ~  e,L  and  use  the  Residue  Theorem  to  integrate 


(2) 

(3) 


pj  +  i<jj  = 


(l-h2-  k2)3!2 
2k 


J-,  (1  +  hsi 


.ijL 


sin  L  +  k  cos  L )2 


dL 


31 


(4) 


Pi  =  (i  +jB)(-byc](k,k) 

*i  =  (1  +jB)(-byS,(k,h) 

where  b  and  B  are  given  by  (2. 1.4-4)  and  (2.1.6-lb),  and 

cJ{k,h)  +  is}(k,h)  =  (k  +  ihy  (5) 

are  obtained  from  the  recursion  formulas 

CJ+1{k,  h)  =  kCj(k,  h)  -  hSj(k,  h) ,  C0  =  1 

Sj+i{k,h)  =  hCj(k,h)  +  kSj(k,h)  ,  50  =  0  ^ 

The  quantities  (L  —  \)m  may  be  calculated  from  the  expansion 

(L  -  A)m  =  +  E(Km  cos  jL  +  ip^smjL)  (7) 

j=i 


Here  the  first-order  coefficients  k\  and  \j){  are  obtained  immediately  from  (2): 


k\  =  0 
Ki  _  2a: 

1  ~  J~ 

=  -&■ 
3 


(8) 


Higher-order  coefficients  K]m  and  are  obtained  using  the  following  recursion  formulas, 
obtained  by  multiplying  the  series  on  the  right  sides  of  (2)  and  (7): 

*m+l  =  E  r(^Km  -  P*^m) 
k-\  K 


^m+l 


tfUi  = 


7*1  +  E  JTT(ai^Km  - 

3  k=\3  t  K 

+ 1  r<***£‘  -  +  (1  -  «;.)  E 

fc=i  K  k=\3  ~  K 

~~Km  -  E  TTT^+fcV’m  +  P;+**m) 

•7  fc=l  J  f  K 

J-l 


(9) 


+  E  r(<T^+fe  +  +  (1  -  fy)  E  7^-rK-fctfm  -  Pi-**m) 

*=i  *  k=\3  ~  K 

Note  that  if  m  is  odd,  then  ( L  —  A)m  is  antisymmetric  about  perigee.  Therefore  we  have 


^  jjL  -  =  0 

i.  e.  ,  only  even  powers  of  ( L  —  A)  have  nonzero  constant  terms. 


(10) 
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Again  we  may  use  (2.4-10)  and  (2.5. 1-2, 3, 4, 7, 8)  to  obtain  the  short-periodic  variations 
rji.  To  do  this,  it  will  be  necessary  to  convert  the  integrals  over  A  into  modified  Fourier  series 
expansions  in  L.  We  begin  by  supposing  /(A)  has  a  modified  Fourier  series  expansion  in  L 
with  known  coefficients  and  it  averages  to  zero: 


M  J 

/(A)  =  C°  +  £  Vm(L  -  A)m  +  £(CJ  cos  jL  4-  sin  jL) 


m—1 


;=i 


<  /(A)  >=  0 


Using 


jmj\  =  j 

the  following  consequences  of  (2. 1.4-2, 5, 6)  and  (2) 

d\  1  /r \2  (1  -  h2  -  k2f!2  ,  .r  .  ... 

=  =====  =  ( - )  =  ,  2  u  •  r  I 'l - nr  =  cos jL  +  ajSinjL),  (14) 

oL  \/l  —  h2  —  k2\a/  (1  +  ft  sin  L  +  k  cos  Ly  ^ 

the  integral 

(L  -  A)"*+1 


(ID 

(12) 

(13) 


J(L-  \)mdX  =  -iil— L-  +  J(L  -  A)mdL  , 


(15) 


the  identities  (2.5.2-7)  and  expansions  (7),  we  can  convert  the  right  side  of  (13)  into  a  Fourier 
series  expansion  in  L.  Note  in  particular  that  equations  (3)  and  (12)  imply  that  the  constant 
term  C°  in  (11)  must  be  related  to  the  Fourier  coefficients  C',S'  and  V'  by 


M  J 

=  -  £  P"*”  -BcVJ  +5i<r>) 

m=l  7=1 

Note  also  that  we  need  the  following  formulas  for  the  product  of  two  Fourier  series: 
E(^J  cos  jL  +  S]  sin  jL)  ^2(pk  cos  kL  -I-  aksin  kL )  =  |  E{^i r(j)(C3pj  +  S3c tj) 

j=1  k= 1  1  7=1 

+  1(1 E  {V-tpt-Si-lot)  +  l\-\?)^C*lpl  +  S*lol) 

^=1 
J 


(16) 


+ y^(cePi+t + st<rj+t) 


e=i 


cos  jL 


(17) 


+ 


(1  -  Sn)  E  (&-**,  +  S’-tpt)  +  E (-C*‘oe  +  Sj+ePe) 

r=max(7— J,l)  t=  1 

J 


+  ^2{C<crj+t  ~  SlPj+t) 


t=l 


S\QJ 


iL} 
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Here  2 *{j)  is  the  inclusion  operator  defined  by 


™  =  {  l  otLfJi*  '  <18> 

Higher-order  integrals  can  be  computed  using  recursion  formulas  obtained  from  equation 
(2.5. 2-9).  We  summarize  the  conversion  for  the  general  multiple  integral  with  the  following 
notation: 

/  f(X)d\k  =  t/fc°(C<,S<,Z>c)  +  £  W?{C<,S<,V<){L-  A)m 

oc  m=1  (19) 

+  E[t/*(C^<^)cosii  +  V?(Cc,Sc,2><)sin  jL\ 
j=i 

(^/(A)«‘)  =  0  (20) 

Here  the  arguments  (C^,5^,2^)  denote  that  the  functions  £/*,  V* ,  and  W™  depend  upon  the 
coefficients  and  Pm  appearing  in  (11)  through  the  relations 

for  j  >  1: 

i  m  j- 1 

i/;  =  -7  +  +  E  (C'-v» + s'-VO 

1  _  m=l  fcs=m«x(j— J,l) 

+i/-'o)  E(-cj*v, + s’+v*)  +  f;(c‘<7,+l  -  5V,+») 

*=i  t=i 

v;  =  i  liw  +  £  +  (i  -  fj,)  £  (C'-%  -  5y-V*) 

3  "  m=l  fc=max(j- J,1 ) 

+ri'-'(»  E(c,+V*  +  +  £(cv,+*  +  sv,+l) 

fc=l  fc=l 

w;  =  -c° 

pm-1 

wr  = -  (2i) 


for  k  >  1: 


V”  =  -  £  w-'f  *2,  -  Y.i'Jlh  + 

m=l  j= 1 


i  Af+fc-1  J-l 

vu,  =  -7  v;+  e  wrp„+p-6n)'£(ut-,v!+vr'p,) 

J  _  m=l  r=l 

+ + V?+V<  -  HV>+<) 

r=i 
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i  M+k-i  j-1 

v{+i  =  '  vi  +  e  +  (i  -  wE^rV,  -  vrv,) 

J  m=l  /=1 


+  E(^+'/>/  +  tfo+<  +  +  v[o,tl) 

/=  l 


w*+1  =  -U? 


wr+l  = 


^r1 


With  this  notation,  we  can  write  the  mean  element  rates  Ai  and  the  short-periodic 
variations  rji  corresponding  to  (1): 

a ,  =  c?  +  Er>r<  +  Efe+^)  (22) 

m= 1  j=l 

A/ +AT  "f-2  oo 

=  C,°  +  £  DT(L-\)m +Y,(Ci  cos  jL  + Si  s\njL)  (23) 

m=l  ji=l 

where 


M+if+2  oo 

cf  =  -  E  A”<-E«^  +  ^) 

m=l  j=l 


mk '  dtk  '  dtk 


—sj-uuc*  S<  2>‘l  +  Y(k  +  ntdlii/i 

~2a  ,6ln  2(C'’*1’  l,  +  ^(  +  '  »*+>  1+2  (  ft*  ’  dtk  ’  ft*  /) 

«  _  Iv>,c<  s<  T)C, .  f  (zlliv*  f^cf  a^f  a*x>f\ 

s,  n*iK,Si,®i)  +  L  „»+!  v*+i  ^  ft*  ■  ft*  -  a<*  j 


_±,  \lv,(c<  c(  p{' ,  f ,,  I  M 

2a  6  U  ^  )5,,Pl)  +  f^k  +  *'  „*+i  Vfc+2  [  dtk  '  dik  '  dtk  Jj 

-^Ixy+Hm^wr^sirf) 

1  lM+k+2(m)(k  I  (dkCi  dk$i 

In  the  absence  of  explicit  time-dependence,  equations  (23)-(24)  can  be  simplified.  The 
short-periodic  variations  become: 


Vi  =  Cf  +  £  D?{L  -  A)m  +  BC?  cos  jL  +  S{  sin  jL) 

m=  1  j=l 
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where 


C? 

Cl 

si 

D) 

D] 

D™ 


M+2 


=  -  £  0"<£-£(C?ft  +  sto) 


m=l 


=  i£/'(Cf,  Sf,  Pf)  -  j Si,  V‘) 
=  ±V>(C<,  Si,  Vi)  -  j^WCi,  Si,  Vi) 

Ac°+±6Msi,vi) 


(26) 


n 

-hv'-i-n6^ 

— Lp— i - - - 

mn  *  2  am(m  —  l)n 


m  —  2 


2.5.4  General  77,  Expansions  in  X,d 


where 


If  one  or  more  perturbations  are  double- averaged,  then  the  functions  Gt  can  be  written  as  a 
double  Fourier  series  in  the  mean  longitude  A  and  the  perturbing-body  phase  angle  0: 

Gi(a,  h,  k,p,  q,  A,  6,  t)  =  ^[Clm(a,  h,  k,p,  q,t)  cos(jA -m0) +Slm(a,  h,  k,p,q,t)  sin(jA -md)] 

J.m 

(1) 

c"  =  hU’:G'ixie 

C{m  =  ]\  f\  G'  cos  O' A  -  mO)d\d0  (2) 

<S/m  =  /  /  GiS\n(jX  —  m6)dAd0 

lit1  J—ir  J—r 

Using  (2.4-11)  and  (2. 5. 1-2),  we  can  obtain  the  mean  element  rates  At\ 


Ai  =  [ C( m  cos(jA  —  mO)  +  5/m  sin(jA  —  mO)] 

(j,m)€e 


(3) 


We  can  then  obtain  the  short-periodic  variations  ??,  by  integrating  the  Equations  of  Averaging 
(2.5. 1-1)  in  a  manner  similar  to  the  single-averaged  case  in  Section  2.5.1.  Assuming  that  the 
coefficients  Cl m,  SI m  and  the  rotation  rate  6  do  not  explicitly  depend  upon  time,  we  obtain 


r)i  —  l^/m  cos(jA  —  md)  +  S{m  sin(jA  —  m#)] 


(4) 


where 


cr  = 

sr  = 


i 


_  3  71  6j 6 


jn  —  md  [  '  2  a  jn  —  md 

1 


jn  —  md  [ 


Qjm  _j_  3  n  8js  ^.jm 

1  2  a  jn  —  md  1 


(5) 
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2.5.5  First-Order  r]ia  for  Conservative  Perturbations 


For  conservative  perturbations,  it  may  be  advantageous  to  use  an  alternate  solution  for  the 
first-order  short-periodic  variations  r),a  which  avoids  having  to  obtain  the  osculating  rate 
functions  Fia.  If  a  perturbation  a  is  conservative,  then  the  osculating  rate  functions  Fia  can 
be  expressed  as  (2.2-6): 


Fta  = 


6 


J=1 


daj 


(1) 


where  (ai  •••a6)  are  the  equinoctial  elements  (a,  h,fc,p,g,  A),  the  quantities  (a,-, a,-)  are  the 
Poisson  brackets  given  by  (2. 1.8-2, 3),  and  K  is  the  osculating  disturbing  function. 

If  the  perturbation  a  is  single-averaged,  we  can  define  the  mean  disturbing  function  U : 


1  tT 

U  =<  V.  >=  —  /  7£(a,  h,  k,p,  q,  A,  t)d\ 

27T  j  —r 


(2) 


Averaging  both  sides  of  (1),  we  obtain  the  first-order  mean  element  rates  due  to  this  pertur¬ 
bation: 

4  4^,  ,du 

A,a  =  -  2^\Oi,aj)—  (3) 

j=i  aai 

Note  that  the  sum  in  (3)  only  goes  up  to  j  =  5  because  =  0. 

The  first-order  short-periodic  variations  caused  by  this  perturbation  can  be  obtained  from 
a  potential- like  function  S  called  the  short-periodic  generating  function: 


S  =  J{H-  U)d\ 


(4) 


<  S  >  =  0  (5) 

Using  (2.3-26),  (2.4-10),  and  (1 )— (4),  we  can  obtain  the  first-order  short-periodic  kernels  £,a 
by  performing  the  integrals  (2.5. 1-2,  3)  (with  the  subscript  i  replaced  by  ia ): 


f  \dS 
S ia  /  Qj)  n 


From  (5)  it  is  clear  that: 

(U)  =  0 

Combining  (2.1.8-2,3)  and  (6)-(7),  we  obtain: 

t  _  2  dS 
^la  “  n2a  8X 


(6) 

(7) 

(8) 


(9) 
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Substituting  (6)  and  (9)  into  (2.5. 1-5),  we  obtain  the  following  expression  for  the  first-order 
short-periodic  variations  T)ia  in  the  absence  of  explicit  time-dependence: 


Via  — " 


U=i 


dS_ 

da} 


+ 


(10) 


In  the  presence  of  explicit  time-dependence,  the  first-order  short- periodic  variations  j/io  are 
given  by: 


Via 


,  £(-l)‘  A,  ,  3  /3‘S  ..  / 

S“[S(°”“,)3^-/‘*r  +(°”  ’A-. 


3  c 

*1 - J^t6 

ncr 


5*5 

dtk 


d\ 


k-l 


(11) 


2.5.6  Second-Order  77,0/3  for  Two  Perturbations  Expanded  in  A 

In  this  section  we  suppose  that  the  osculating  rate  functions  Fia  and  the  first-order  short- 
periodic  variations  r)ia  can  be  written  as  Fourier  series  in  the  mean  longitude  A: 

Fia(a,h,k,p,q,\)  =  C°a{a,h,k,p,q,t) 

0° 

+  £[C(«»  h,  k,P,  9*  *)  cos  +  C(a>  h,  k,P,  V,  0  sin  ;A]  (1) 

j=i 

Via  =  ^2(Cja  cos;  A  +  5^  sin; A)  (2) 

i=i 

From  (2.3-27),  the  second-order  functions  GIO/?  are 


/-i  _  v'  ^Fia  _  ,  15n  f  _  dvi, 

G.o/J  =  ~^TTlrP  +  fo2<>i6VlaVl0  ~  2L 


r=l  d°r 


lT0 


(3) 


Substituting  (1)— (2)  into  (3),  and  using  (2.5.3-17)  with  «/  =  00,  we  can  write  Giap  as  a 
Fourier  series  in  A: 

Gia0  =  C/3  +  £(C/3  cos;  A  +  S/a/3  sin  ;  A)  (4) 

j=i 
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where 


CL , 


=  £ 
r=l 


dC?  ;  0  dC\a 

r0~  r0lh7 


1  -Si 


daT 


■  J~1  / 
aE( 


~d^~Cr0~ 


dsj;k  k 
~d^~br0 


) 


-jQ0sia  + 


1  -6 


ii 


i- 1 

E 

^=1 


(i  -  *)($ 


•i-fc 


r1* 

u6/3 


+  ci 


s“> + 


'IP 


Then  the  second-order  mean  element  rates  are 


Aiap  —  Ciap 

and  the  second-order  short- periodic  variations  are 


Via0  —  Y,(Cfa0 cos j\  -(-  S3ia0sinjX) 
j= i 


(6) 

(7) 
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r 


where  the  C3a0  and  S3ia0  are  given  by  (2.5.1-15)  (with  the  subscript  i  replaced  by  ia(i)  in 
terms  of  C\a0  and  S3a0.  The  formulas  in  this  section  were  published  in  [Danielson,  March 
1993]. 

2.5.7  Second-Order  T)ia0  for  Two  Perturbations  Expanded  in  L 

In  this  section  we  suppose  that  the  osculating  rate  functions  Fia  and  the  first-order  short- 
periodic  variations  Tjia  can  be  written  as  finite  modified  Fourier  series  in  the  true 
longitude  L : 

Fia(ai  h,k,p,q,  L,t)  =  Cfa(a,h,k,p,q,t) 

Ja  (1) 
+  £l£/a(a> h'  P.  9)  0  cosJL  +  SL(aih,k,p,  q,  t)  sinjl] 

i=i 

Ma 

Via(at  h,  k,p,q,  L,i)  =  Cia(a,  h,  k,p,q,  t)  +  ^  Dia(a,h,k,p,q,t)(L  —  A) 

Ka  m=1  (2) 

+DO».m.p.*o  cos  kL  -f  Ska(a,  h,  k,p ,  q,t )  sin  kL] 

k-l 

The  second-order  functions  Gia0  are  again  given  by  (2.5.6-3).  Differentiating  (l)-(2),  we 
can  obtain  the  needed  partials 


Jar 


=  C°:  +  52(C3;  cos  jL  +  S3:  Sin  jL) 

dn  AfQr  KaT 

+  £  DTQr(L  -  A)m  +  £«£  cos  kL  +  S£  sin  kL) 

U°r  m=l  fc=l 


(3) 


(4) 


The  product  of  two  Fourier  series  can  be  converted  into  a  single  Fourier  series  with  the 
formula 

J  K 

Y,(C]  cos  jL  +  S3  sin  jL)  53(C'fc  cos  +  Sk  sin  kL) 

j=i  fc=i 

i  J+A' 

=  o  E  {zT'{J'K)(j)(cjcj  +  s3s3) 

1  j= i 

mm(j-l,K)  min 

+[ij+K(j)  £  (c3~kck -s}-ksk)+it~l(j)  53  {cj+kck  +  sj+ksk) 

k=m&x(]-J,\ )  k=l 

min(A— j,J) 

+1 ?~\j)  53  ( CkCj+k  +  SkSj+k)\  cos  jL 

*=1 

min(j-X,A)  min(J— j,K) 

+[t/+*0)  53  (c3-ksk+s3~kck)+i(-1(j)  53  (_Ci+^  +  5i+*c*) 

k=aact(j-J,l)  k=l 

min(A'  -  j,J) 

+Z,K-'U)  E  (C*S'+*-S‘C>+*)]  sin  ji}  (5) 

*=1 
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With  the  use  of  (2. 5. 3-7),  we  then  obtain 

Mar+Mff  oo 

Gia0  =  C°a0  +  £  V?a0(L-  A)m  +  £(<&,  cos  jL  +  S}%a0  sin  jL)  (6) 

m=l  }=l 

where 


.  min {JaT,K*) 

+;  E  (CtC‘a  +  srste)  -  A,eC” 
i  >=■ 


cu 


Z{Z,K‘U) czcie + iroxscj, 

r=I 

+  lir+^0)  £  «£*'<%  -  $iak'TSr0) 

fc=max(j-Jar,l) 

,  min  (Jar-j,Kff) 

+k r-'u)  z  (c{,+*’'cr*9+5i+‘^) 

/  jt=i 

,  min(A^-j,yor) 

+5ir-‘(»  E  «t<?4+‘  +  s£sS‘) 
z  *= 1 

M* 

-ir'U)AeC£  +  E  flSPTO'XS*!!. 

m=l 

+i(i  -  <;.)  E  «£■*■’*:,  - 

Z  k=max(j-J°r,l) 

+\ir~'U)  Je'(C+‘-'^  +  5L+‘”*)  +  5  E«£<+‘  +  5‘r«,+‘)]} 

z  *=1  z  *=I 
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+ i ru)c?0cu 


So2 

j  _  vain(j-lJ<|,) 


(ci;kci,  -  si;kst.) 


i-k  ck 


+^rkK°u)  e 

Z  k=roAx(j—Ka,l) 

i  min  (Ka-jj<0) 

+5xr-‘(j)  E  (c£‘c,‘»  +  s;:‘sf„) 

z  fc=l 

,  min  (K*-j,K°) 

+2  E  (<?,*„<?;♦*  +  s‘.s;f) 

L  k=  1 

M*  i  j- 1 

+  E  fl&PT y)c?.4  +  ;0  -  *>.)  E  (cfc***  -  sfcV* ) 

m=l  4  fc=maji(j-/fa,l) 

+kr-o)  eV;«+**‘  +  s£V£) + i  E(c*,<+‘ + s,‘0«,+*)] 

*=i  *  *=i 


+  E^[xrO')cia«”+i(i-«;,)  E  (c;;‘4  -  s;jv* ) 

m=l  Z  i=nui|;-A’0,l) 

Kp-} 


E 

+kK”-o)  eW*‘  +  sftVi) + ;  £;«??„<+* + sf„«+*)]} 


*=i 


(7) 


Stop 


EWuvzs*,  +  jru)s;:c^e 

T—\ 

(S>J-k,T  qk  ,  qj  —  k,T  y-lk  \ 
(Via  JTp  +  °ia  Is rp ) 


i  min(j— 1  ,Ke) 

+^r*K,u)  e 


*=maxO-Jar,l) 

,  min(Jar— 

+5xr-(i)  e  (-c+‘-'s;fl+5i+*-'c^) 

z  fc=i 

1  min  [Kfi-j,Jar) 

+5  xf'-‘c>)  E  «tsjT‘  -  5£c;;*) 

z  fc=i 

M* 

-lTU)A,esr  +  £  Br„(xr'(j)5'>“ 

m=l 

^  *=m*xO-Jar,l) 

+ jxr-'o)  JeVci;‘>‘  +s'„+‘"4) + J  E(C«,+l  -  s£<+‘))} 


*=1 


k=\ 
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+ zro)o; 


8a2 


la 


+|  ir*K’u)  E  (<&% + s; :‘cf„) 

fc=m*x(j— K"0,!) 
nun  (Ka-j,Ke) 


+kr-l(j)  z . (-«:% + s;:‘c;s) 


4=1 

••  min  (K*—j,Ka) 

+kr-'(j)  E  (OiJ‘  -  &<#*) 


k= 1 


Af'’  i  j-i 

+  e  +  5<i  -  E  (ccvi + «;;*£.) 

m=l  k=max(j-Ka,l) 

+\ir-,U)Z\-c;tk^  +  s£‘*:L)  +  5  E^L’tit1'  -  5,‘„<+‘)] 

Z  4=1  Z  4=1 

+  EDr.[zfO)s^  +  5(i-«i.)  E  +  s;j‘<4) 

m— 1  k=max(j- K0,l) 

+\z?‘-'U)  e'(-c;JV‘  +  +  i  E(cs,^‘  -  sk<+‘)]} 


4=1 


4=i 


r=l 

+S^[;r"”(m)C?'*0W  +  Zf'»Cj,I>r.  +zr+*"(m)  E 

00  J=1 


m— 1 


Then  the  second-order  mean  element  rates  are 


Mar+M0  oo 

■4«=c?„a+  E  *>;>“+ E(0>i + s; >j) 

m=l  j=l 


(8) 


and  the  second-order  short-periodic  variations  are 

Af^+M^+K+a  oo 

*«>  =  C?„„  +  E  -  *)m  +  Z(c’.emiL  +  SL„ sin jL)  (9) 

m=X  j=l 

where  the  Cja0  ,  S\a0  ,  D%0  are  given  by  (2.5.3-24)  (with  the  subscript  i  replaced  by  ia/3 
and  M  replaced  by  Mar  +  M0)  in  terms  of  C\a0  and  S-a0.  The  formulas  in  this  section  were 
published  in  [Danielson,  August  1993]. 


2.6  Partial  Derivatives  for  State  Estimation 

Observational  data  may  be  used  to  improve  the  estimate  of  a  satellite’s  state.  Some  differ¬ 
ential  correction  algorithms  which  have  been  used  in  conjunction  with  SST  are  described  in 
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[Green,  1979],  [Taylor,  1982]  and  [Long,  Capellari,  Velez,  and  Fuchs,  1989].  It  is  our  purpose 
here  only  to  explain  how  to  obtain  the  partial  derivatives  needed  in  such  filters. 

We  let  Ok  denote  the  value  of  the  kth  observed  quantity  computed  with  the  SST  orbital 
generator.  Tie  SST  state  variables  are  the  initial  mean  elements  a,{t0)  and  various  constant 
parameters  c,  (the  geopotential  coefficients  Cnm  and  5„m  in  (2.7-1),  the  drag  coefficient  Co 
in  (3.4-3),  the  solar  radiation  pressure  coefficient  Cr  in  (3.5-6),  etc.  ).  Required  for  a  batch 
filter  are  the  partial  derivatives  of  the  Ok  with  respect  to  the  state  variables  a, (to)  and  ct. 

The  actual  observations  are  commonly  of  position  and  velocity  components  in  a  local 
coordinate  frame  fixed  on  the  surface  of  the  Earth.  However,  through  transformations  these 
components  may  be  expressed  in  terms  of  the  orbital  elements,  so  we  can  regard  the  Ok  as  an 
implicit  function  of  the  osculating  elements  a.,.  Application  of  the  chain  rule  then  produces 


dOk  A  dOk  dhj 
dai(t 0)  do*  dai(to) 

(1) 

dOk  A  dOk  da, 
dci  da.)  dcx 

(2) 

Assuming  we  can  obtain  the  partials  analytically,  our  remaining  task  is 

partials  and  |^.  Differentiating  the  decomposition  (1-1)  yields 

to  calculate  the 

dtj  _  y*  (c  ,  drb  \  dak 
da{{t 0)  Aj  \  lk  dak)  da, {t0) 

(3) 

Mb  _  y*  (s  4-  dak  i  dr]} 

dc,  k=i\jk  dak )  d^  dc. 

(4) 

The  partials  and  are  often  arranged  to  form  a  matrix  $,  referred  to  as  the 

state  transition  matrix: 


*  dai 

dai(t0) 


*(Mo)  = 


das 

dai(to) 

0 


0 


dai 

dai 

dai  ' 

das{to) 

dci 

dct 

das 

da6 

das 

da6{to) 

dci 

dc( 

0 

1  0 

...  o 

J 

0  1 

‘ 

; 

l 

1  0 

0 

0  ••• 

0  1 

(5) 


Here  £  is  the  number  of  parameters  c*.  Differentiating  (5)  with  respect  to  t  and  interchanging 
the  ordinary  and  partial  derivatives,  we  can  obtain  the  following  initial  value  problem  for 
*(Mo): 

*(Mo)  =  F*(Mo),  *(Wo)  =  I  (6) 
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Here  I  is  the  identity  matrix  and 


F  = 


dd\ 

dai 

dd\ 

das 

0 

ddi 

dc\ 

dd\ 

dct 

ddo 

da6 

0 

dd6 

dd6 

dax 

das 

dcx 

dct 

(7) 

0 

•  0 

0 

... 

. 

o 

The  $  matrix  is  a  function  only  of  the  five  slowly  varying  mean  elements,  and  therefore 
the  numerical  integration  of  (6)  can  be  done  with  the  same  large  step  size  as  used  in  the 
integration  of  equations  (l)-(2)  for  the  mean  element  rates.  Values  of  $  at  observation  times 
not  coinciding  with  the  integrator  step  times  can  be  obtained  by  interpolation. 

Our  task  has  thus  been  reduced  to  the  calculation  of  the  partial  derivatives  |^,  |^, 

and  jh..  These  same  partials  are  also  needed  in  a  sequential  Kalman  filter.  For  the  two-body 
part  ol  the  mean  element  rates 


the  only  nonzero  partial  in  (7)  is 


Oi  —  TlSfO 

(8) 

dde  3  n 

dax  2  a 

(9) 

and  thus 


$  = 


1 

0 

0 

0 

0 

3 n(t  -  t0) 
2a 
0 


0 

1 


0  1 


(10) 


0 


1  0 
0  1 


Although  the  differential  correction  algorithm  for  updating  the  initial  mean  elements  may 
converge  with  only  the  two-body  partials  (10),  the  speed  of  convergence  can  be  improved 
by  including  the  dominant  perturbation  partials.  Analytical  formulas  have  been  obtained 
for  the  partial  derivatives  with  respect  to  the  mean  elements  a,  of  the  J2  contribution  to 
the  mean  element  rates  dj  (equations  (3.1-12)),  the  partial  derivatives  with  respect  to  the 
geopotential  coefficients  Cnm  and  Snm  of  the  resonant  tesseral  contribution  to  the  mean 
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element  rates  a :  (equations  (3.2-9)),  and  the  partial  derivatives  with  respect  to  the  mean 
elements  a,  of  the  J2  contribution  to  the  short  periodic  variations  rj}  (from  the  expressions 
in  Section  4.1). 


2.7  Central-Body  Gravitational  Potential 

The  well-known  expression  for  the  disturbing  function  due  to  the  gravitational  field  of  the 
centred  body  is  [Battin,  1987]: 


N  min (n,M)  /  Dv  n 

=  -  J3  53  (-)  Pnm(sin <(>)(Cnm  cos mrp  +  5„m  sin m^)  (1) 

r  n=2  m=0  ' r  ' 


Here 


r 


<t> 

0 


R 


Cnm,  Sn 


M 

N 


radial  distance  from  center  of  mass  of  central  body 

geocentric  latitude 

geographic  longitude 

central-body  gravitational  constant 

central-body  mean  equatorial  radius 

associated  Legendre  function  of  order  m  and  degree  n 

geopotential  constant  coefficients 

maximum  order  of  geopotential  field  (M  <  N) 

maximum  degree  of  geopotential  field 


In  this  section  all  elements  are  osculating  (even  though  they  do  not  have  hats). 

In  the  first  subsection  we  shall  outline  the  development  of  the  central- body  gravitational 
potential  into  the  form  used  in  SST.  Complete  details  are  to  be  found  in  [Cefola,  1976], 
[McClain,  1978],  [Proulx,  McClain,  Early,  and  Cefola,  1981],  and  [Proulx,  1982].  Then  in 
the  remaining  subsections  we  describe  methods  for  calculating  the  various  functions  used  in 
the  expansion. 


2.7.1  Expansion  of  the  Geopotential  in  Equinoctial  Variables 
We  start  by  writing  (2.7-1)  in  the  complex  form 


{  N  mi n(nM)  ,  p>.n 

-E  E  (-)  em(.in^)(c. 

r  n=2  m=0  ^  r  ' 


lSnrn)c 


•mV' 


(1) 


Here  i  =  y/—\  and  Re  {z}  is  the  real  part  of  z.  With  the  goal  of  expressing  (1)  in  terms  of 
the  equinoctial  elements,  we  set 

^  =  aB  -  9  (2) 


where  0  is  the  central  body  rotation  angle  and  a g  is  the  right  ascension.  If  we  let  (xg,  yB,  z b) 
denote  a  right-handed  orthonormal  triad  fixed  in  the  central  body,  with  Xb  pointing  to  the 
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prime  meridian  and  Zg  to  the  geographic  north  pole,  then  0  may  be  calculated  from  [Early, 
1982]: 


sin  0 


cosO 


-f  •  ya  +  /g  •  xg 
1  +  /7 

f  •  xg  +  /g  •  yB 
1  -f  I'f 


(3) 

(4) 


(Remember  that  7  is  defined  by  (2.1.9-lc).) 

Next  the  spherical  harmonics  P„m(sin  <t>)e'maB  are  expanded  as  a  Fourier  series  in  the 
true  longitude  L,  using  a  rotational  transformation  theorem  for  the  spherical  harmonics: 


F„m(sin <j>)eim°B  =  £  V?aS™t''L 

j=— n 

The  Vn”)  coefficients  are  defined  by: 


ym 
v  ns 


f  (-l)"2’  (n  +  s)\{n  -  s)! 

2"  (n  -  m)!(af*)!(a=*)! 


if  n  —  s  is  even 


0 


if  n  —  s  is  odd 


(5) 


(6) 


The  rotation  functions  S™3(a,f3,i)  may  be  expressed  in  terms  of  the  dot  products  (a,/?, 7) 
of  the  zg  vector  with  the  equinoctial  reference  triad  (f,g,  w): 


(_  l)m-*2*(a  +  i/3)lm~3{  1  + 

(-1)m-.2-m^  ~V”w"  ~  "!f  (o + upr-'v + h 

(n  -f  s)!(n  —  s)! 

2_1(a  -  1  +  h)Imp:zT’s+m(i) 


(7) 


if  .s  <  — m 
if  |s|  <  m 

if  s>  m 


m 

(Remember  that  /  is  defined  by  (2. 1.2-2).)  Here  Pevw( 7)  are  Jacobi  polynomials.  (Note  that 
commas  are  used  to  separate  indices  in  a  symbol  such  as  PZ+s*'~m~*  in  order  to  prevent 
ambiguities.) 

Next  the  product  (1 )ne'aL  is  expanded  in  a  Fourier  series  in  the  mean  longitude  A  (the 
sixth  equinoctial  element  a6): 


e«jA 


Now  the  Hansen  coefficients  X™  are  defined  by  [Hansen,  1855] 


OO 

,,f  =  Y,  xTe 

j=-oo 


ijM 


(S) 


(9) 
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and  the  kernel  A-"*  of  the  Hansen  coefficient  X™  is  defined  by 


K*\e)  =  e~is-’'X™(e)  (10) 

where  /  is  the  true  anomaly,  M  is  the  mean  anomaly,  and  in  (10)  e  is  the  orbital  eccentricity. 
Hence,  remembering  equations  (2. 1.2-1,  4),  we  can  express  the  functions  YJna  as 


Y^(h,  k)  =  [k  +  ih  sgn  (s  -  i)]|4-j|  A"4 


(11) 


Here 

1  if  x  >  0 
—1  if  x  <  0 

The  last  step  is  simply  a  rearrangement  of  the  order  of  summation,  so  as  to  isolate  the 
total  disturbing  potential  due  to  the  phase  angle  j A  —  m0 ,  and  to  facilitate  the  use  of  stable 
recursion  formulas.  We  also  introduce  some  new  notation,  to  enable  the  results  to  be  written 
concisely.  First,  define  the  functions 


f  (— l)m-424(l  +  I')) 


-Im 


if  s  <  —m 


C,(7)  =  (-11 - 2— <?  + / "ili"  -  "ftl  +  h)“  if  W<m 

(n  +  s)!(n  —  s)! 

2-4(l  +  I'i)Im  if  s>m 


Next,  put 

G]m .  +  iHi 


-{ 


[A:  +  ih  sgn(s  —  j)]!*  ;*(q  +  iJ0)m  h  if  |s|  <  m 


[A:  +  ih  sgn(s  —  j)]l*  ^[a  —  ifi  sgn(s  —  m)]l*-/ml  if  |s|  >  m 

Then  define  the  Jacobi  polynomial  PfW  indices  by 

—  m  if  |s|  <  m 
|s|  if  |s|  >  m 

v  =  \m  —  s| 
w  =  |m  +  s| 

The  disturbing  function  can  now  finally  be  written  as 


l 


(13) 


(14) 


(15) 


M 


N 


N 


,R, 


EE  E  Qnimv™rzKr~l'apr 

j=-oo  m= 0  a=  —  N  n=max(2,m,|s|) 

(Gin.  + 


(16) 


Note  that  the  functions  G^s  and  defined  by  (14)  are  of  degree  |s  —  j  |  in  the  eccen¬ 
tricity;  the  power  |s  —  j\  is  called  a  D’Alembert  characteristic.  In  SST  the  number  of  terms 
retained  in  the  expansion  (16)  is  truncated  by  prescribing  the  maximum  possible  value  of 
the  D’Alembert  characteristics,  together  with  the  maximum  degree  N  (and,  if  desired,  the 
maximum  order  M  <  N)  of  the  geopotential  field. 
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2.7.2  Calculation  of  Coefficients 

The  V™  coefficients  are  defined  by  (2.7. 1-6).  Since 


i/m 

n,s 


(-i)JKT_3 


(i) 


we  can  restrict  our  discussion  to  the  case  s  >  0  without  loss  of  generality.  Furthermore, 
since  Vn"  =  0  when  n  —  s  is  odd,  we  need  only  consider  the  case  when  n  —  s  is  even.  Also, 
note  that  the  lowest  value  of  the  degree  n  in  the  summations  (2.7.1-16)  is  greater  than  or 
equal  to  2,m,  and  )s). 

Suitable  recurrence  relations  are 

(n  -r  s  -f  l)(n  -  s  +  1) 

(n  -  m  +  2)(n  -  m  +  1)  (2) 

=  (n-m)VZ 

Appropriate  initialization  is  provided  by 

1 

(2^  +  1)  Q  (3) 

(a  +  1) 

I.e.,  to  calculate  the  coefficients,  first  use  (3)  to  get  values  for  m  =  0  and  n  =  s  for 
s  =  0, 1,  •  •  •  .  Then  use  (2b)  for  m  >  0  and  still  n  =  s.  Finally,  use  (2a)  for  increasing  n 
with  any  nonnegative  m  and  s. 


V° 

vo,o 


V° 


VI 


fl+2,3 


T/m+ 1 

r  n,3 


2.7.3  Calculation  of  Kernels  Kfa  of  Hansen  Coefficients 

(Throughout  this  section  e  represents  the  eccentricity.)  From  the  definitions  (2. 7. 1-9,  10), 
the  kernels  of  the  Hansen  coefficients  are  given  by 

[w  cos(5/  - jM)dM  (o 

Z7T  J-1 r  a 


The  kernels  of  the  Hansen  coefficients  are  thus  functions  of  the  orbital  eccentricity  e.  Note 
from  (1)  that 


Kn,.  =  Knr 


n,—  8 


(2) 


so  we  can  restrict  our  discussion  to  the  case  s  >  0  without  loss  of  generality. 

For  the  special  case  j  =  0,  the  kernels  may  be  evaluated  in  a  form  algebraically  closed 
in  the  eccentricity.  This  is  because  the  Hansen  coefficients  with  j  =  0  are  related  to  the 
associated  Legendre  functions  by  [McClain,  1978] 


vn»5 

^0 


(n  -  l)!xn 

(n  +  s  -  1)! 
(-l^n-s  +  l)!*-"-1 


(n+1)! 


pn-i(x)  for  n  >  1 

(x)  for  n  >  0 


pi 

rn+W 


(3) 
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where 


vr^  vi  -  w  -  ^  b  v  7 

and  many  recursion  formulas  are  available  for  the  associated  Legendre  polynomials,  which 
for  arguments  in  the  range  1  <  x  <  00  are  defined  by 

1  jn+» 

K(x)  =  ^(x2  -  i)'/2^(x2  -  D”  (5) 

For  the  special  kernels  with  the  first  superscript  negative  (needed  in  Sections  3.1  and  4.1), 
appropriate  recursion  formulas  are 


AT"'1’'  = 


for  n  =  s  >  0 
for  n  =  s  +  1  >  1 


f(2n  _  3) I<on'a  -  (n  -  2)Kon+1'3]  for  n  >  s  +  2  >  2 
(n  +  5-  l)(n  —  s  —  1)  1  J 


For  the  special  kernels  with  the  first  superscript  nonnegative  (needed  in  Section  3.2),  appro¬ 
priate  recursion  formulas  are 


K'*  = 


with  initializations 


— —Ao" 2,3-1  for  n  =  s  -  1  >  1 

s 

{2s  +  1ll(Z~u  for  n  =  s  >  1 

s  +  1  0 

2n  +  1  ryn_ha  (n  +  s)(n  -  s)  r^n_2,s  r__  _  ^  ^  0 

— A0  -  .  „  K0  for  n  >  s  -f  1  >  2 

n-fl  n(n  -f  l)xz 


for  n  =  s  >  1 


Kq°  =  1 

8) 

K1  =  -1 

The  general  kernels  A-"-1’3  may  be  computed  from  the  following  recurrence  relation 
[Proulx,  McClain,  Early,  and  Cefola,  1981]: 

Kr'- = (3--n)(1-483)(i-n-,)  {(3  - n)(1  -  ")(3  -  w 


■(2  -  n)[(3  -  n)(l  -  n)  +  y  ]  A';”+,'‘ +  j!(l  -  n)A7"+i’j 


To  initialize  this  recurrence  relation,  we  need  the  values  of  the  three  kernels  Kj  n’3  ,  K}  n+1’3, 
and  K~n+3'a  at  n  =  max(2,  m,s).  These  latter  kernels  are  calculated  by  infinite  series 
representations.  In  a  study  of  the  various  possibilities  [Proulx  and  McClain,  1988],  it  was 
found  that  the  expansion  of  choice  is 


K?  =  (1  -  e!)"+}  £ 
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Here 


(11) 


a  =  ma x(j  —  s,  0) 
b  =  max(s  —  j,  0) 

The  Y™  terms  are  called  modified  Newcomb  operators  and  may  be  computed  by  the  recur¬ 
rence  relation 


4(/<  +  <’)C‘  =  2(2»  -  n)V2£*  +  (>  -  n)Y?Z‘  -  2(2S  +  n)Y^^ 


n,j+2 


^n,j—  1 


~(s  +  n)Ypy:22  +  2(2 p  +  2a  -  4  -  n)Y?l\<<r_x 


(12) 


Recursion  formula  (12)  is  initialized  by 

y0y  =  i  (i3) 

and  by  treating  quantities  with  negative  subscripts  as  identically  zero.  Note  that  the  Y™ 
terms  are  rational  constants,  and  therefore  they  can  be  computed  once  and  stored  for  all 
later  applications. 


2.7.4  Calculation  of  Jacobi  Polynomials  P™ 

The  Jacobi  polynomials  appear  in  the  expression  (2.7. 1-7)  for  the  rotation  functions.  PfW( 7) 
is  a  polynomial  of  degree  £  in  7,  which  from  (2.1.9-lc)  is  the  cosine  of  the  angle  between 
a  vector  from  the  geographic  south  to  north  pole  of  the  central  body  and  the  angular  mo¬ 
mentum  vector  of  the  satellite.  The  Jacobi  polynomials  P?w( 7)  with  m  =  0  in  the  indices 
(2.7.1-15),  i.e.  £  —  n  —  s  >  0  and  v  =  w  =  s  >  0,  are  related  to  the  associated  Legendre 
functions  Pn„( 7)  by 

(1) 

The  Jacobi  polynomials  can  be  computed  from  the  standard  recurrence  relation  [Szego, 
1959]: 

2£{£  +  v  +  w)(2£  +  v  +  w-  2)P(W(~f)  = 

(2£  +  v  +  w-  1)[(2£  +  v  +  w){2£  +  w  +  to  -  2)7  +  v2  -  w 2}P?”M  (2) 

-2{£  +  v-  l)(£  +  w-  1){2£  +  v  +  w)P??2{ 7) 

This  recursion  formula  is  initialized  by 


(3) 


2.7.5  Calculation  of  G3m>  and  H3ma  Polynomials 

From  the  definitions  (2.7.1-14),  the  functions  G3ma  and  H 3ma  are  polynomials  in  the  equinoctial 
elements  h, k  and  the  direction  cosines  a,/?. 
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These  polynomials  may  all  be  calculated  from  one  set  of  generic  recurrence  formulas, 
based  on  the  Cj  and  Sj  polynomials  obtained  from  (2.5. 3-6): 


G>mt 


HL 


C\s-j\(k,  h)Cm-i3(a,fl)  -  Isgn(s  —  j)S^ji(k,  h)Sm.ja{a,/3)  for  |s|  <  m 

C\s-j\{k,h)C\a-im[(a,p)  +  sgn(s  -  j)sgn{s  -  m)S\s-:\(k,h)S[3-lm\(a,p)  for  \s\  >m 

(1) 

h)Sm-Ia(a,  0)  +  sgn(s  -  j)S\,-j\(k,  /3)  for  |s|  <  m 

-sgn(s  -  m)C\a-j\(k,  h)S|,_/m|(<*,  $)  +  sgn(s  -  j)S\a-j\(k,  h)C|,_/m|(<*,  0)  for  \s\  >  m 

(2) 


2.8  Third-Body  Gravitational  Potential 


The  disturbing  function  due  to  the  gravitational  field  of  a  third-body  point  mass  is  [Battin, 
1987]: 


7l(r,  <f>,  t) 


H3  (  R3  rcos(j)\ 

R3\\K.3-r\~-lir) 


(1) 


Here 

r  =  vector  from  the  center  of  mass  of  the  central  body  to  the  satellite 
R3(<)  =  vector  from  the  center  of  mass  of  the  central  body  to  the  third  body 
<f>  —  angle  between  the  vectors  r  and  R3 
H3  =  third-body  gravitational  constant 

In  this  section  all  elements  are  osculating  (even  though  they  do  not  have  hats). 
The  quantity  can  be  expanded  in  the  following  series: 


R3 

|R3  -  r| 


1 


2  r  cos  <t>  1  r2 

R3  +w3 


S  ikY  K(cos*] 


(2) 


where  Pn  is  the  Legendre  polynomial  of  degree  n.  Hence  the  third-body  disturbing  function 
(1)  can  be  written  as 

R=l|(i)"p"{cos«  (3) 

where  N  is  the  maximum  power  of  the  parallax  factor  ^  to  be  retained  in  the  expansion. 

The  further  development  of  the  third- body  potential  into  the  form  used  in  SST  is  similar 
to  that  of  the  central-body  potential  in  Section  2.7.  Complete  details  are  to  be  found 
in  [Cefola  and  Broucke,  1975],  [McClain,  1978],  [Cefola  and  McClain,  1978],  and  [Slutsky, 
1983]. 


2.8.1  Expansion  of  Third-Body  Potential  in  Equinoctial  Variables 

With  the  goal  of  expressing  (2.8-3)  in  terms  of  the  equinoctial  elements,  we  set 

cos^  =  acosL  +  flsmL  (1) 
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where  L  is  the  true  longitude,  and  now  (a,/?,  7)  are  the  dot  products  of  the  unit  vector 
^  with  the  equinoctial  reference  triad  (f ,  g,  w).  Due  to  the  motion  of  the  third  body, 
(a,  /?,  7)  are  slowly  varying  functions  of  the  time  t  and  are  the  source  of  weak  time-dependence 
effects. 

Next  the  expression  Pn(a  cos  L  +  Ps'm  L)  is  expanded  into  a  Fourier  series,  using  an 
addition  formula  for  the  Legendre  polynomials: 

Pn(acosL  -f-  jdsinL)  =  J^(2  ~  S0t)Vn,Qnt(^)[C3(a,  P)  cos  sL  +  S,(a,fl)  sin  sL]  (2) 

4=0 


Here  60s  is  the  Kronecker  delta,  and  the  polynomials  Ct(a,0)  and  5,(o, /?)  are  the  same  as 
those  in  (2. 7.5-1,  2).  The  new  coefficients  introduced  in  (2)  are  defined  by: 


V«4  = 


(~1  )V  (n~s) ! 

m?  (*?)'■ 


2n 

0 


if  n  —  s  is  even 
if  n  —  s  is  odd 


(3) 


Q»4<7)  = 


d‘Pn(  7) 

dY 


After  substituting  (2)  into  (2.8-3),  we  can  write  the  result  in  the  complex  form 

^  =  Re  ED  2  -  So.)  (jj)“  Vn.Qn.[C.(cl,fi)  -  iS,(a,0 )]  Q”  e"L| 


(4) 

(5) 


The  last  steps  are  simply  the  replacement  of  the  product  elsL  by  the  expansion 
(2.7. 1-8),  and  a  rearrangement  of  the  order  of  summation.  The  disturbing  function  can  then 
be  written  as 

R  =  R<- 1  £  £  £  £  (2-«o.)(-|-lVn,y'’“l3n.[C.(a,/?)-t5.(a,^)]ei'4  (6) 

l  U3  i=-<»4=0n=max(2,4)  J 


2.8.2  Calculation  of  Vn3  Coefficients 

The  Vn„  coefficients  are  defined  by  (2.8. 1-3).  Since  Vnt  =  0  when  n  —  s  is  odd,  we  need  only 
consider  the  case  when  n  —  s  is  even.  (Note  from  (2. 7.1-6)  that  V™  =  Vna .) 

A  suitable  recurrence  relation  is  [Cefola  and  Broucke,  1975] 


Ki+2  ,«  — 


(n  -  s  -f  1) 
(n  +  s  +  2) 


Appropriate  initialization  is  provided  by 


K),o  =  1 

K+ 1,4+1  = 


(1) 

(2) 
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2.8.3  Calculation  of  Qna  Polynomials 

The  Qna  polynomials  defined  by  (2.8. 1-4)  are  derivatives  of  the  Legendre  polynomials  eval¬ 
uated  at  7,  which  from  (2.2-7c)  is  the  cosine  of  the  angle  between  R3  and  the  angular 
momentum  vector  of  the  satellite.  The  polynomial  Qna  can  also  be  expressed  in  terms  of  the 
associated  Legendre  function  Pna: 

QM  =  (i  -  (i) 

(Note  from  (2.7.4-1)  that  />“. (7)  =  VffoQM  ■) 

Recursion  formulas  for  the  Qng  polynomials  follow  directly  from  standard  recursion  for¬ 
mulas  for  the  Pna  functions  [Cefola  and  Broucke,  1975]: 


QnM 


(2s  -  l)Q4_lii_1(7) 

(2s  +  lb<?«,,(7) 

(2n  -  l)7<?n-i,a(7)  -(n  +  s-  l)Qn.2,ah) 
(n-s) 


for  n  =  s 
for  n  =  i  +  l 

for  n  >  s  • f  1 


These  recursion  formulas  are  simply  initialized  by 


(2) 


Qo,o  —  i 


(3) 


3  First-Order  Mean  Element  Rates 

As  we  have  seen,  the  first-order  mean  element  rates  A,c  are  given  by  equations  (2.4-18).  The 
osculating  rate  functions  F1Q  for  a  conservative  perturb ~cion  are  given  by  (2.2-10)  and  for  a 
non  conservative  perturbation  by  (2.2-5).  In  this  section  we  record  the  specific  form  of  these 
equations  for  each  of  several  perturbations. 

3.1  Central-Body  Gravitational  Zonal  Harmonics 

For  the  central-body  gravitational  zonal  harmonics,  the  appropriate  averaging  operator 
<  •••  >  is  (2.4-10)  and  the  appropriate  disturbing  function  Tt  is  (2.7.1-16)  with  m  =  0. 
Further  details  of  the  reduction  of  the  averaged  equations  to  the  form  recorded  here  may  be 
found  in  [Cefola  and  Broucke,  1975]  and  [McClain,  1978]. 

The  first-order  contribution  of  the  central- body  gravitational  zonal  harmonics  to  the 
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averaged  equations  of  motion  (2.4-1)  is 


da 

dl=° 

dh  BdU  k  ,  rr  . 

dt  “  Adk  +  ABipi’ay  IqU'P'l) 
dk  _  BdU  h 
dt~  Adh  AB{pU,cn  IqU,(h) 


dp 

dt 

<h 

dt 


2AB 

1C 


Ui(h 


(1) 


2  ABUun 


dX  _  2 adU_ 

dt  A  da  ^ 


B  /LdU  ,dU.i  1  ,  „  , 

A(  1  +  B){h  dh  +  k  dk  ]  +  AB{pU,ay  IqU'01  ] 


Here  (a,h,k,p,q,  A)  are  now  the  mean  elements  and  U  is  the  mean  disturbing  function 
(2. 5. 5-2).  In  deriving  (1)  from  (2.2-10),  we  have  made  use  of  the  following  property  of  the 
cross-derivatives  for  the  mean  disturbing  function: 


U fhk  U,Q(3  —  0 


The  mean  disturbing  function  reduces  to 


u  ^  JL  / R\n 

=  E(2-«o.)(-)  JnVn.Ko"-'’-Q„,G. 


Here 


Jn  —  CnO 


Kna 

oLw 

G.=  G°o, 


geopotential  coefficients 

coefficients  calculated  from  (2.8.2- 1,2) 

kernels  of  Hansen  coefficients  calculated  from  (2. 7.3-6) 

functions  calculated  from  (2. 8. 3-2, 3) 

polynomials  calculated  in  Section  2.7.5 


Since 

G.  +  iH.  =  G? ,  +  ./&  =  (*  +  ih)‘(a  -  iff)- 


an  alternate  set  of  ursion  formulas  for  the  G„  polynomials  is 


G,  =  (ka  +  h0)G,.x  -  (ha  -  k0)Ha.x  ,  G0  =  1 
Ha  =  (ha  -  k0)G,.x  +  (ka  +  hff)H„-X  ,  H0  =  0 


(2) 

(3) 


(4) 

(5) 


Note  that  the  Ga  polynomials  are  of  degree  s  in  the  eccentricity;  for  small  eccentricity  orbits, 
the  series  (3)  may  be  truncated  by  prescribing  the  maximum  possible  value  of  the  D’Alembert 
characteristic  s. 
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In  equations  (1),  we  need  the  partial  derivatives  of  U  with  respect  to  (a,  h,  k,  0,^,7). 
These  are  easily  obtained  by  differentiating  equation  (3): 


dU_ 

da 

dU 

dh 

dU 

dk 

dU 

da 

dU 

dp 

dU 

d'f 


4  E(2  -  M("  +  l)(-rJ«Vn,K;"-''Qn,G, 

“  s,n  ** 

--  E(2  -  <b.)(-)V.  v„<3„,  +  <>X3G 

a  V  °h 


dKo 


-n-l,s 


dx 

dU  H  c  Y  Q  ( T/-n-l,tdG»  .  ,  ldK0  \ 

dU  dG, 


) 


~  =  -^E(2-M(-)vnv;^0-n-^Q 

/I  T  AT 


3,n 


«,n 


a 

=  -^£(2  -W(-)V. v„k;"-'"Q 

"  a 


dG, 

!  ns  rv 

oa 

dG, 
'  d/3 


(6) 


--  E(2  -  M(-)”4K,A-0-”-’'‘^c, 

^  "  a  d~f 


3,n 


Here  we  have  obtained  the  partial  derivatives  with  respect  to  h  and  k  of  K0  "  1,3  (x)  from 
the  chain  rule  and  the  definition  (2. 7. 3-4). 

Recursion  formulas  for  — fl~</x  ~  are  obtained  by  differentiating  the  recursion  formulas 
(2. 7. 3-5): 


dKon~hs 


dx 


=  < 


(l+2s)X 


2s 


(«  -  i)x2 


(n  +  s  —  l)(n  -  s  -f  1) 
+  2  n_M 

X 


j  is~n<a 

(2n  —  3)— y - (n  —  2) 


for  n  —  s 
for  n  =  s  +  1 
dKon+h3 


dx 


(7) 


dx 

for  n  >  s  +  1 


Recursion  formulas  for  are  obtained  by  differentiating  (2. 8. 1-4): 


dQn,> 

d'i 


(7)  =  Qn,,+l(~t) 


(8) 


Recursion  formulas  for  the  partial  derivatives  of  G,  may  be  obtained  by  differentiating  (4): 


dG, 

dh 

dG, 

dk 

dG, 

da 

dG, 

dd 


s/lG,-!  -  saH,.  1 
sqG,_i  -f  s/3H,-i 
skG,- 1  —  shH,- 1 
shG,- 1  +  skH,- 1 


(9) 
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If  we  retain  only  the  J2  term  in  the  expansion  (3),  the  central-body  gravitational  dis¬ 
turbing  function  simplifies  to 

(.0) 


U  = 


where 


a3(l  -h1-  fr2)3/2 
.  3  pR2J2 


The  contribution  of  (10)  to  the  averaged  equation  of  motion  (1)  is 


dh  _  Jfc[372  -  1  -I-  27 (pa  —  Iq/3)\ 

dt  AB4a3 

dk  _  Jh[ 372  —  1  -f  27 (pa  —  Iq/3)] 

It  =  AB4a? 

dp  CJfa 

Tt  =  -  AB4  a3 

dq  ICJa'i 

dt  AB4a3 

dX  _  J[(  1  +  B)(372  -  1)  +  27  (pa  -  IgP)\ 
dt  AB4a3 


In  order  to  update  the  orbital  elements  in  a  differential  corrections  procedure,  it  is  nec¬ 
essary  to  compute  the  partial  derivatives  with  respect  to  the  mean  orbital  elements  of  the 
mean  element  rates  (the  |^-  in  the  matrix  F  defined  by  (2.6-7)).  The  partial  derivatives  of 
the  J2  contribution  (11)  to  the  mean  element  rates  are  easily  obtained,  using  (2. 1.6-1)  and 
(2. 1.9-4).  The  nonzero  derivatives  are 
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dh 

da 

d'h 

dh 

d'h 

dk 

dh 

dp 

dh 

dq 

dk 

da 

die 

dh 

dk 

dk 

dk 

dp 

die 

dq 


7Jk[3-y2  -  1  +  2 lipc*  -  Iq0)\ 

2a4AB4 

4Jhk[3~f2  -1-1-  27 (pa  -  Iq/3)] 
a3AB 6 

J(l  -  h2  +  3fc2)[372  -  1  4-  2i{pa  -  IM 
a3AB 6 

2Jfe[6a7  +  2 p(a2  -  72)  -  2q2a'f  -  2 Iqfi(a  +  p7)  +  Cq7] 

a3AB4C 

21Jk[6/3 7  -f  2pa7  +  2/97*  -  2Ipqct'i  -  2p2/?7  +  C07) 

a3AB4C 

7Jh[3~j2  -  1  +  27 (pa  -  /g/?)] 

2a4  AB4 

(1  -  fc2  +  3fe2)J[372  -  1  +  2^ (pa  -  /gg)] 
a3AB« 

4Jhk[3-y 2  —  1  +  2/y{pot  —  /g/?)] 

a3A5« 

2  J/t[6a7  4-  2p(a2  -  72)  -  2g2a7  -  2IqP(a  A  pi)  4-  Ca7] 

a3AB4C 

2/Jft[6/?7  4-  2pa7  +  2/^72  -  2 Ipqa-j  -  2p2/?7  4-  C/37] 

a3AB4C 
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dp 

ICJfa 

da 

2a4  AB4 

dp 

\CJhfa 

dh 

a3AB6 

dp 

iCJkfa 

dk 

a3AB* 

dp 

2J\p0'r  +  <*(/?  +  /?7)1 

dp 

a3AB4 

dp 

2IJ[-ICqfa  +  /32  -  72  +  po7] 

dq 

a3AB4C 

dq 

ICIJa-i 

da 

2a4  AB4 

dq 

\CIJha-) 

dh 

a3AB 6 

dq_ 

ACIJka 7 

dk 

a3AB6 

dq 

21J\pa~j  +  a2  —  72  —  Iqflj] 

dp 

a3AB4 

dq 

2  JlqCIa'y  +  P(a  -  pj)] 

dq 

aMB'C 

d\ 

7J[(l  +  B)(37J-l)+2l(pa-/9«] 

da 

2  a4AB4 

dX 

JM(372  -  1)(4  +  5 B)  +  8 7(pa  -  Iq0)] 

dh 

a3AB 6 

dx 

JA:[(372  -  1)(4  +  35)  +  87(pa  -  /<?/?)] 

dk 

a3AB6 

dX 

2J[6(1  +  B)a~f  +  2  p(a2  —  72)  —  2  q2a7  —  2Ipq(3~i  +  Ca7] 

dp 

a3AB4C 

dX 

2/J[6(l  +  B)&~i  +  2p(3(a  -  P7)  +  2/97(7  -  P«)  +  CV?7] 

dq 

(12) 


59 


3.2  Third-Body  Gravitational  Potential 

For  a  third-body  point  mass,  the  appropriate  averaging  operator  is  (2.4-10)  and  the  appro¬ 
priate  disturbing  function  71  is  (2.8. 1-6).  Further  details  of  the  reduction  of  the  averaged 
equations  to  the  form  recorded  here  may  be  found  in  [Cefola  and  Broucke,  1975]  and  [Mc¬ 
Clain,  1978]. 

The  first-order  contribution  of  the  third-body  gravitational  disturbing  function  to  the 
averaged  equations  of  motion  is  identical  in  form  to  equations  (3.1-1)  for  the  central-body 
zonal  harmonics.  Of  course,  the  direction  cosines  (a,  /?,  y)  have  different  interpretations  for 
the  two  perturbations. 

The  mean  disturbing  function  is  now 


(1) 


Here 

Vn3  =  coefficients  calculated  from  (2.8.2- 1,  2) 

Kqs  =  kernels  of  the  Hansen  coefficients  calculated  from  (2. 7.3-7,  8) 

Qns(l)  =  polynomials  calculated  from  (2. 8.3-2,  3) 

Ga  =  polynomials  which  can  be  calculated  from  (3.1-5) 


The  partial  derivatives  of  U  needed  in  equations  (3.1-1)  are  easily  obtained  by  differentiating 
equation  (1): 


% = t  p  ■ -  «*>  (£)"  K-Q-  (K?w + h^-i¥) 

f = f  »2  -  *•>  (£)"  Ww + k^G-iJP) 

-u^ur^. 


aG, 

da 

dGs 


(2) 


dKn~ 

Recursion  formulas  for  ■  are  obtained  by  differentiating  the  recursion  formulas 
(2. 7.3-6,  7): 


dK y 

<*X 


=  < 


0 


n— 1,» 


2n  +  1  dK% 
l  n  +  1  d\ 


for  n  =  s  —  1  or  n  =  s 

m—2,3 


(n  +  s){n  -  s)  dKp~  2(n  +  s)(n  -  s)  2,s 
n(n  +  l)x2  d\  n(n  +  l)x3  ° 
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n(n  +  !)x3 


for  n  >  s 

(3) 


Recursion  formulas  for  the  derivatives  of  Qnt  and  G,  are  given  by  (3.1-8,  9). 

If  we  wish  to  retain  only  the  dominant  terms  in  the  expansion  (1),  we  should  include  at 
least  the  first  two  terms,  since  the  n  =  2  term  vanishes  in  the  ^  and  ^  equations  for  zero 
eccentricity  orbits,  leaving  the  n  =  3  term  dominant  [Collins  and  Cefola,  1978]. 


3.3  Central-Body  Gravitational  Resonant  Tesserals 

For  the  central-body  gravitational  tesseral  harmonics,  the  appropriate  averaging  operator  is 
(2.4-11)  and  the  appropriate  disturbing  function  is  (2.7.1-16)  with  m  ^  0.  Further  details  of 
the  reduction  of  the  averaged  equations  to  the  form  here  may  be  found  in  [Proulx,  McClain, 
Early,  and  Cefola,  1981]  and  [Proulx,  1982]. 

The  first-order  contribution  of  the  central-body  gravitational  tesseral  harmonics  to  the 
averaged  equations  of  motion  (2.4-1)  is  identical  to  (2.2-10),  except  that  77  is  replaced  by  U: 


a 

h 

k 

P 

9 

A 


hB  dU 


A(l  +  B)  d\\ 


2 adU_ 

A  d\ 

BdU  k  (  \ 

\BdU  h  (  rr  T  rr  \  kB  dt/1 
.Adh  +  AB\pU'a'r  IqU^)  + 

C 

2AB 
C 

2  AB 

2a  dU  B 
~  A  da  +  A(  1  4-  B) 


p(v,u  u~-,)  d\) 

i- 

i 

(1) 

j 

-  W,a, 

f.du  ,dU\  1  /  IT  r  tr  \ 
\hdh+kdk)  +  AB  IqU ^  ) 


Here  (a,  h,  k ,  p,  q ,  A)  are  now  the  mean  elements  and  U  is  the  mean  disturbing  function  defined 
by 


U  =  <  ^  >  =  /  /  ^(ai  h,  k,  p,  q,  A,  0,  t)dX  d9 

+Re  I  2^2  E  £  /'  R(«,  M,  J>,  <,,  A',  9',  ,)e-iO'A <(«'] 

(2) 

Remember  that  B  here  is  the  set  of  all  ordered  pairs  (j,  m)  with  the  properties  (2.4-12,  13). 

The  mean  disturbing  function  is  identical  to  (2.7.1-16),  except  that  now  only  the  resonant 
tesserals  are  included  in  the  summations: 


,  n  m  N  N 

t^MjEE  E  E 

j  m=l  s——N  n=max(2,m,|j|) 

(j,Tn)€B 


(*)m  i-vzTvcr^pr 

(GL  +  iHL)(Cnm  -  iSnm)e^-^j 


(3) 
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In  equations  (1),  we  need  the  partial  derivatives  of  U  with  respect  to  (a,  h,Jb,  A,  0,^,7). 
These  are  easily  obtained  by  differentiating  equation  (3): 

fr  =  R'{-jr  E  <"  +  0  (7)" i'"vzrZK;'-'-pr(Oi.  + 


QTJ  a  /  R\n 

~  =  R*{;  E  (7)  rvzTZPr(c„-iS„) 


(aai,  mi 

[K>  ("aT  +  'Sir 

%  -  M7  E  (f)"  imKK.f7'(c„n  -  is„m) 


+  2h(G‘m,  + 


9GL.  1  .9HL 

dk  8k 


dK~n~ l’4 

+  2k(Gi,  +  H’m.)  ' 


—  =  Re{^  ■£,  3  (7)  C'V”C,A--"-1-'fT(Gi„  +  iH’m.)(Cm  -  (4) 

§7  =  M£  E  (f)"''mcr:,ft-r",''fr(^i+'^!i)(c„„-is„„)ei«i-",i} 

~  =  Re(;  E  (^)Vv”i7.A7*',’,Pr(^i  +  i^i)(C--i5™)e,(ii-”s|} 

%  =  Re{“  E  (7)’ ImVZKr-l*(G’m,  +  iHU{C„m-iSnm) 

j*wi,5,n 

(Arm  APVW\  x 

Here  we  have  obtained  the  partial  derivatives  with  respect  to  h  and  k  of  /^^n_1,'(e2)  from 

the  chain  rule  and  the  relation  h2  -f  A:2  =  e2,  where  e  is  the  orbital  eccentricity. 

dKnt  J 

Recursion  formulas  for  are  obtained  by  differentiating  the  expansion  (2.7.3-10): 

dKns  (n  4-  31  oo 

sk = -|rr$A7' + 0  -  *2>"+l  E  w 

Recursion  formulas  for  are  obtained  by  differentiating  the  recursion  formulas 

(2. 7.4-2,  3): 
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2  i(£  4  v  4  w)(2i  +  v  +  u>- 2)—~(j)  = 

«7 


jpvw 

\art-\ 


(2(  4  v  4  w  —  l)[(2f?  4  v  4  w)(2(  4 ■  v  +  tv  —  2)7  4  v2  —  w 2] — 7—^(7) 

07 

dPvw 

—2(1  +  t;  -  1)(*  4  w  -  l)(2i  4  u  4  w)—j=l( 7) 

07 

4(2^  4vrii'  -  1)(2^  4  v  4  w)(2^  4  u  4  to-  2)P™1(i) 

dPvw 

^-=0 

d'l 


Recursion  formulas  for  the  partial  derivatives  of  GJm3  and  H3mi  are  obtained  by  differentiating 
(2.5.3-51  and  (2. 7.5-1,  2): 

\s  -  j\C\a-j\-i(k,h)Cm-u{ot,0)  -  l(s  -  j)S\a-j\-X(k,h)Sm-is(a,p)  for  |s|  <  m 
=  \s  -  j\C\s-j\-i(k,h)C\s-im\[a,P)  -  (s-j)sgn(m  -  s)5|*_i|_1(fc,  h)S\3.Im\(a,  P) 

for  |s|  >  m 

(7) 

etc.  Formulas  for  are  obtained  by  differentiating  (2.7.1-13). 

In  order  to  update  the  geopotential  coefficients  Cnm  and  Snm  in  a  differential  corrections 
procedure,  it  is  necessary  to  commute  the  partial  derivatives  with  respect  to  the  coefficients 
of  the  mean  element  rates  (the  in  the  matrix  F  defined  by  (2.6-7)).  These  are  easily 
obtained  by  partial  differentiating  (1)  and  (4)  with  respect  to  Cnm  and  Snm.  Introducing 
the  parameter 

fl  if  max(2,  m,  |s|)  <  n 

’  (0  otherwise  '  ' 

we  can  write  the  results  in  the  compact  form 


da  .  da  2 pi 

dCnm  +  'ds„„  A 


E  E  o(^y  rv?j;,K-”-'‘pr(G!n,  +  iH’m,y'’ 

j=-oo  5=- N  VG/ 


If  we  assume  that  the  orbital  eccentricity  is  zero,  the  averaged  equations  of  motion  for 
the  resonant  tesserals  significantly  simplify.  See  [Collins  and  Cefola,  1978]. 


3.4  Atmospheric  Drag 

For  atmospheric  drag,  the  appropriate  averaging  operator  is  (2.4-10),  and  the  first-order 
mean  element  rates  are  obtained  by  substituting  (2.2-5)  into  (2.4-18).  To  avoid  having  to 
solve  Kepler’s  equation,  and  to  smooth  the  perturbation  around  perigee,  we  can  convert  the 
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integrals  over  the  mean  longitude  A  into  integrals  over  the  eccentric  longitude  F  or  true 
longitude  L  by  use  of  (2.5. 2-6)  or  (2.5.3-14).  The  first-order  contribution  of  drag  to  the 
averaged  equations  of  motion  can  then  be  written  in  either  of  the  forms 


da 

It 


or 


dai 

~dt 


1  /Ls  fr\2  (da,  \  jj 

-  h2  -  k2  JLi  \aj  \  dr  ) 


(la) 

(16) 


The  quantities  ^  and  ^  are  given  in  terms  of  the  equinoctial  elements  by  equations 
(2. 1.4-1,  6,  7,  8,  9),  (2. 1.6-1),  and  (2.1.7-3). 

The  limits  (F\,  F2)  in  (la)  or  (Li,L2)  in  (lb)  indicate  the  values  of  F  or  L  at  atmosphere 
entry  and  exit.  If  the  satellite  enters  and  leaves  the  atmosphere  at  a  critical  distance  r  from 
the  center  of  the  central  body,  then 


F\  =  —  E  u>  IQ 
F2  =  E  uj  ICt 


where 


or 


where 


E  =  arccos 


L\  =  — ■/ 

L2  =  f  +  u>  +  IQ 


f  =  arccos 


(26) 


Of  course,  if  the  satellite  remains  totally  within  the  atmosphere,  the  limits  of  integration  in 
(1)  can  be  taken  to  be  (— 7r,7r). 

The  perturbing  acceleration  due  to  atmospheric  drag  is  commonly  modeled  by  the  formula 
[Escobal,  1965]: 

q= vl(r- v)  (3) 

Here 

Co  =  drag  coefficient  of  satellite  ( Cd  =  2  for  a  sphere, 

Cd  =  4  for  a  flat  plate  perpendicular  to  v) 

A  =  cross  sectional  area  of  satellite 
m  =  mass  of  satellite 
p  =  density  of  atmosphere 
f  =  ^  =  velocity  of  satellite 
v  =  velocity  of  atmosphere 

If  we  assume  that  the  atmosphere  rotates  with  an  angular  rate  equal  to  the  angular  velocity 
uj  of  the  central  body,  then  v  =  w  x  r.  The  vector  q  is  resolved  along  the  (x,y,2)  axes  of 
Figures  1  and  2  for  use  in  the  quadratures  (1). 


The  developers  of  SST  have  used  various  density  models  for  the  upper  atmosphere  of  the 
Earth.  One  of  these  is  the  modified  Harris-Priester  atmosphere  (described  in 
[Long,  Capellari,  Velez,  and  Fuchs,  1989]): 

P  =  Pmin  d"  (Pm*x  Pram)  (t- )  (4) 

where 


Pmin  (H) 
Pmax(H) 
Hmin 

Hrnixx 


Hj -H 

Pmim  e 

Hi -H 

Pm&x\  e  »">»• 

H2-H2 

lnf52^) 

yPminj  J 

HX~H2 

In  (^2.) 

\Pm*xj  / 


Here  pmin(H )  and  pm^iH)  denote  the  minimum  and  maximum  densities  at  a  height  H  above 
a  reference  ellipsoid  ( Ht  <  H  <  H2),  Hi ,  H2,  pm im,  Pm**,,  Pmin2?  Pmax2,  N  are  constants 
whose  values  are  available  from  Tables,  and  fa  denotes  the  angle  between  the  diurnal  bulge 
and  the  satellite.  If  b  denotes  a  unit  vector  pointing  from  the  Earth’s  center  to  the  diurnal 
bulge,  then 

b  •  r 

cos  fa  = -  (5) 

r 

The  diurnal  bulge  follows  the  path  of  the  Sun  but,  because  of  the  Earth’s  rotation,  lags  the 
sub-solar  point  by  an  angle  6b  (approximately  30°  at  2  P.M.  local  time).  We  can  obtain  the 
vector  b  by  rotating  the  vector  Z#  (pointing  from  the  Earth’s  center  to  the  sun)  through  an 
angle  6b  about  the  Earth’s  axis  of  rotation.  Letting  R  denote  the  3x3  matrix  whose  elements 
are  direction  cosines  between  the  (a:,  y,z)  axes  and  an  Earth-fixed  set  of  Cartesian  axes,  we 
can  write  the  transformation  law  between  the  (x,y,z)  components  of  b  and  z b  as 


(6) 


«x' 

COS0fc 

—  sin  6b 

0  ' 

'  ZBx  ' 

by 

=  R 

sin  6b 

COS  &b 

0 

Rr 

*By 

K  J 

0 

0 

1 

Z Bz 

Here  Rr  denotes  the  transpose  of  the  matrix  R  (see  [Danielson,  1991]). 


3.5  Solar  Radiation  Pressure 

The  general  equations  expressing  the  first-order  contribution  of  solar  radiation  pressure  to  the 
averaged  equations  of  motion  are  formally  identical  to  the  equations  (3.4-1)  for  atmospheric 
drag. 

The  limits  (Fx,F2)  in  (la)  or  (Li,L2)  in  (lb)  now  indicate  the  values  of  F  or  L  at  shadow 
exit  and  entry.  If  we  assume  the  shadow  is  a  circular  cylinder  in  shape,  the  shadow  exit  and 
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entry  longitudes  are  determined  by  the  shadow  equation  (as  explained  in  [Escobal,  1965]  and 
[Cefola,  Long,  and  Holloway,  1974]): 

5  =  0  (1) 

Here 

5  =  1  —  M{  1  +  kcosL  +  hs\nL)2  —  ( acosL  +  flsin  L)2 


M 

a 

P 


Rl 


a2(l  -  h?  -  A:2) 

r3 

TT  g 


R: 


where  is  the  central-body  radius,  and  R3  is  the  vector  from  the  center  of  mass  of  the 
central  body  to  the  sun.  To  obtain  the  solutions  to  equation  (1),  the  following  quartic 
equation  must  be  solved: 

Ao  cos4  L  +  A\  cos3  L  +  A2  cos2  L  +  A3  cos  L  +  A4  =  0  (2) 


where 


Ao  =  4  B2+C2 
Ax  =  8B  Mh  +  AC  Mk 
A2  =  —AB2  -f  AM2h2  —  2P  C  +  AM2k2 
A3  =  -8B  Mh-AVMk 
A4  =  —AM2h2  +  V2 
B  =  a0  +  Mhk 
C  =  a2  -  02  +  M{k2  -  h2) 

V  =  1  -  P2  -  +  h2) 


The  real  roots  of  (2)  must  be  sorted  to  eliminate  extraneous  roots  and  to  determine  the  exit 
and  entry  values  of  true  longitude  L.  The  correct  values  of  L  must  satisfy  (1)  as  well  as  the 
condition 


R3 

R3 


•  -  =  cos  (f>  =  a  cos  L  +  P  sin  L  <  0 


(3) 


At  exit  from  shadow 

dS 

dL 

while  at  entry  into  shadow 

dS 

dL 


>  0 

(4) 

<0 

(5) 

Of  course,  if  the  satellite  remains  totally  within  sunlight,  the  limits  of  integration  in  (3.4-1) 
can  be  taken  to  be  ( — 7r,  7 r). 
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The  perturbing  acceleration  due  to  solar  radiation  pressure  is  [Cefola,  1982]: 


_  CrA  C  2  (r  -  R3) 
q_  2m  c  ®|r-R3|3 


(6) 


Here 

Cr  —  radiation  pressure  coefficient  of  satellite 

(Cr  =  2  for  a  spherical  mirror  or  black  body, 

Cr  —  4  for  a  flat  mirror  perpendicular  to  (r  —  R3)) 

A  =  cross  sectional  area  of  satellite 
m  =  mass  of  satellite 
C  =  mean  solar  flux  at  surface  of  sun 
c  =  speed  of  light 
Rq  =  radius  of  sun 

r  —  R3  =  position  vector  from  sun  to  satellite 
If  we  suppose  that  the  satellite  is  always  in  sunlight,  and  that  the  parameter 

(7) 

is  constant,  we  can  derive  (6)  from  the  disturbing  function 

(8) 

Use  of  the  expansion  (2.8-2)  then  leads  to 


K  = 


|R3  -  r| 


CrAC  2 
"  2m  c  0 


n=-ki(wJp^  (9> 

The  radiation  disturbing  function  (9)  is  identical  to  the  third-body  disturbing  function 
(2.8-3),  except  that  the  factor  /x3  is  replaced  by  — T  and  the  summation  starts  from  n  =  1. 
Hence  we  can  immediately  write  down  the  mean  radiation  disturbing  function  by  analogy 
with  (3.2-1): 

v  =  ~E  £  (2-«(j)"w»c.  (10) 

3  3=0  n=max(l,s)  ® 

If  we  retain  only  the  first  nonzero  term  in  the  expansion  (10),  the  mean  radiation  dis¬ 
turbing  function  simplifies  to 

U  =  V(ka  +  h/3)  (11) 

where 

3Ta 


The  contribution  of  (11)  to  the  averaged  equations  of  motion  (3.1-1)  is 
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dt 
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( kp  ~  Ihq) 
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~  ^jj(kp-Ihq) 


(12) 


4  First-Order  Short-Periodic  Variations 

Knowing  Fourier  series  expansions  for  the  osculating  rate  functions  Fia  or  the  osculating 
disturbing  function  7 Z,  we  can  construct  the  first-order  short-periodic  variations  rjia  from  the 
results  in  Section  2.5.  In  this  chapter  we  record  the  specific  forms  of  the  expansions  for  each 
of  several  perturbations. 


4.1  Central-Body  Gravitational  Zonal  Harmonics 

For  the  central-body  gravitational  zonal  harmonics,  the  appropriate  disturbing  function  71 
is  (2.7.1-16)  with  m  =  0.  From  the  results  in  Sections  2.5.1  or  2.5.5,  we  can  construct  a 
Fourier  series  expansion  in  the  mean  longitude  A  for  the  first-order  short-periodic  variations 
t)ia,  as  was  done  by  Green  [1979]. 

However,  for  the  central- body  zonal  harmonics,  it  is  possible  to  construct  a  finite  modified 
Fourier  series  in  the  true  longitude  L  for  the  r/1Q.  For  single- averaged  perturbing  forces  which 
increase  rapidly  as  the  satellite  approaches  the  central  body,  notably  central-body  zonal 
harmonics  and  atmospheric  drag,  a  Fourier  series  expansion  in  L  will  converge  faster  than 
equivalent  expansions  in  F  or  A  when  the  eccentricity  of  the  satellite  orbit  is  large.  This 
happens  because  the  magnitude  of  such  a  perturbation  is  strongly  peaked  around  perigee 
when  considered  as  a  function  of  A.  Transforming  the  independent  variable  to  L  broadens 
the  peak  considerably,  making  it  easier  to  approximate  with  a  finite  Fourier  series.  In  this 
Section  we  outline  the  construction  of  this  most  desirable  expansion  in  L.  Further  details  may 
be  found  in  [Cefola  and  McClain,  1978],  [Kaniecki,  1979],  [McClain  and  Slutsky,  1980],  and 
[Slutsky,  1980]. 
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The  disturbing  function  less  it’s  mean  can  be  written  as 


R  -  V  =  Re  {-£  £  £(2  -  So.)  (-)"  J„Vn.Q  « .[C.W)  +  iS.(a,0)\  f;  Y~n~l’~MeijX  1 

(  a  n=2i=0  j-- oo  J 


i*o 


(1) 


Here 


Jn  =  —Cno  =  geopotential  coefficients 

Vna  =  coefficients  calculated  from  (2. 8. 2- 1,2) 

Qna( 7)  =  polynomials  calculated  from  (2. 8. 3-2,  3) 

Ca(a,ft),  Sa(a,0)  =  polynomials  calculated  from  (2. 5. 3-6) 

y-n-i.-s  _  coefficients  0f  the  expansion  (2.7. 1-8): 

(n  \  n+l  00 

-J  e~iaL  =  £  Y~n~1’~lei^x  (2) 

r'  j=- 00 

The  short-periodic  generating  function  (2.5. 5-4,  5)  is  easily  obtained  by  integrating  (1): 

{.,  N  n  /R\n  00  y-n_1'-4e‘J>  ^ 

a  n=2  »=0  V  °  J  j=- 00  lJ  ) 

i#  0 

(3) 

The  infinite  series  in  the  mean  longitude  A  in  (3)  can  be  replaced  by  a  finite  modified 
series  in  the  true  longitude  L.  To  see  this,  first  integrate  both  sides  of  (2)  with  respect  to  A 
to  obtain 

OO  yr-n— 1,— *  ij\  .  /a\n+l 

E  J  ...  -  /  (~)  e-‘^X-XY0 


-n-1  ,-a 


j=-oo 


(4) 


Next  perform  the  integral  in  (4)  by  using  the  expansion 


=  x/l  -  h2  -  k 2  Y0-n~2’-jeij 
and  the  change  of  variable  (from  2.5.3-14) 


(5) 


j=-n 


dx  =  ■  1  —  (-Y 

x/l-h2-  k2  W 


dL 


(6) 


The  infinite  series  then  becomes 


j=-oo 


j=-(n-l) 


(7) 
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Replacing  the  infinite  series  in  (3)  with  (7),  and  introducing  the  kernels  K0  n  lj  of  the 
Hansen  coefficients  through  (2.7.1-11),  we  obtain 


S  =  U(L-  A) 


+  iS.to.flHCmlM)  -  t»gn(j)5m(t,  k)}Kr-lJ^1-’>L 


} 


(8) 


Since  all  of  the  symbols  in  (8)  except  for  i  =  \/—\  are  real,  we  can  easily  cast  S  into  the  real 
form 


S  «  1/(1  -  A)  -  £  £  »2  - 

°  n=2  3=0 


ff  *.-—J 

w=® 

+  nf^on-W 

i= i 


//j,  cos(j  -  s)L  -f  Gj,  sin (j  -  s)L 


J  -s 

-Ija  cos  (j  +  s)Z,  +  Jja  sin(;  +  $)L 


j  +  s 


where 


1) 

(9) 


Ln(i)  =  (-)  v;,gn,(7) 

Gj.  =  C,{k,h)Ca{a,0)  +  Sj{k,h)Sa(a,  0) 

H}a  =  Cj(k,h)Sa(aJ)-Sj(k,h)Ca(a,0) 

Ij3  =  Cj(k,h)St(a,0)  +  SJ(k,h)Ca(a,0) 


(10) 


J]a  =  Cj(k,h)Ca(a,  0)  —  Sj(k,h)Sa(a,0) 

The  last  step  is  a  redefinition  of  indices,  so  as  to  isolate  the  coefficients  of  cos  jL  and 
sin  jL,  and  a  rearrangement  of  the  order  of  summation.  The  short-periodic  generating  func¬ 
tion  can  then  finally  be  written  as  the  finite  modified  Fourier  series 


2N-i 


S  =  C°  +  U{L  -  \)+  £  ( Cj  cos  jL  +  Sj  sin  jL) 

:=i 


Here  the  coefficient  C°  is 


2/V-l 

C°  =  -  ’£(C‘Pi  +  S’ai) 
1=1 


where  pj  and  cr}  are  given  by  (2.5.3-3,  4).  The  coefficients  are 

C>  =l?(j)Cjl  +C* 

where 

Ci'  =  ~^7  )  fs;  E  (2  -  6o,,-i)J.H.,..jK^n-l  ‘L-„- 

aJ  (  L»=i  n=»+i 


(ii) 


(12) 


(13) 
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(14) 


*=0  n=m»*(j+ij+l) 


-  2 


C1 2 


=  a{ 

»j  l 


z?0) 


E(2  -  <0 j-.W..j-.KZ‘-ULr 

3=  1 

E  E  (2  - 

«=ln=;+l 


+  I|'v-10)E!j 


(15) 


j-l  minO-l.A/') 


Ei  =  E  E  (2-h,1-.)Jj.,j-.Ko  ‘Lr 

a=j~  min(j— 1,7V)  n—j—a 
min(j-l,/V)-l  min(j-l,Af) 

+  E  E  (2  -  for  j  even 


(16a) 


»=2 


n=a+l 


min(j— l,iV)  —  1  min(j-l,A/) 

=  E  E  (2-S0,j.a)JJ,,j-,Kon-l’sLr 

*=^i  n=*+1 

(  0  for  AT  =  2,3 

for  j  odd 

Ap  min(j-l,A?) 

2|"-30)  E  E  (2  -  for  w  >  4 

«=j-min(j-  1,/V)  n=j—a 

(m 

Similarly,  the  coefficients  SJ  are 

S]  =  If  {j)S]1  +  Sj2 


+  < 


(17) 


where 


s'1  =  -Alif'-’o)  e  E  (2-v.-j)^g.,-)k0— '■•/.r' 

“7  (  [*=■?  n=4+l 


(18) 

+  E  E  (2  -  foj+.)JnG.j+.Kr'''L?'  +2Z?U)JjGojK;i-,*Li\ 


*=0  n=max(j+*,jH 


5j2  = 


-4  k"(i) 

“j  l 


s=l 

i  iV 


+zT,(;)  E  E  (2 - 6oj-.)JnJ..i-,Ko"-‘',Li- 

3=1 n=j+l 


+  l!iv-,0)Ej| 


(19) 
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(20a) 


n  = 


£  — 1  min(j— X,Af) 

E  E  (2-v.)U,-.C'ir' 

a=j-mia(j-l,N)  n=j-a 
min(j— 1,N)— 1  min(j-  1,JV) 

+  E  E  (2-  for  i 

._i  n=»+X 


min(j— 1,JV)—  1  min(j-l,N) 

25  -  E  E 

t=l=l  n=*+l 

0  for  N  =  2, 3 

for  j  odd 

If-3U)  E  E  (2  -  6oj-.)J«J.j-.KZn-l*Li:*ioT  N  >  4 

,  a—j— min(j— 1,N)  n=j—s 

(206) 

Note  that  the  first  index  of  the  G,  //,  /,  J  polynomials  defined  in  (10)  indicates  their  degree 
in  the  eccentricity;  for  small  eccentricity  orbits,  the  series  (14)— (16)  and  ( 18)— (20)  may  be 
truncated  by  prescribing  the  maximum  possible  value  of  the  D’Alembert  characteristic  s. 

The  first-order  short-periodic  variations  i)ta  generated  by  the  function  S  given  by  (11) 
can  be  derived  using  equations  (2.2-10),  (2.5.2-T),  (2.5.5-10),  and  the  following: 

dL  1,3,,  f[  /  h2  +  k2\  n/n  ,2„ 

^  =  —jp{2kb+kbyl  +  — 2 — J+2  {B +  k2b)  cos  L 

+2hkbsinL  +  —[B  +  (A:2  —  62)6]cos2Z/  +  ^( B  +  2k2 b)  sin2L| 


%  -  UlhB+hb(i+ 


h2  +  k2' 


+  2hkbcos  L 


+2 (B  -f  h2b)  sin  L  +  —  [— B  -1-  (A:2  —  62)6]  cos  2 L  4-  —  {B  +  2h2b)  sin  2 L\ 

Li  &  * 

%r  —  7T~~ - f  2A: cos L  +  2h sin L  +  —  „  cos2L  +  hksin2L\ 

aX  B* 1  2  2  J 

Here,  as  usual,  h  and  k  are  equinoctial  orbital  elements  and  6  and  B  are  defined  by 
(2. 1.4-4)  and  (2.1.6-lb). 

In  the  absence  of  explicit  time-dependence,  the  first-order  short-periodic  variations  can 
be  written  as  the  finite  modified  Fourier  series 


Via  =  Cf  +  Di(L  —  A)  +  ^2  (C?  cos  jL  +  S{  sin  jL)  (22 

j= i 

(Remember  that  the  equation  of  the  center  ( L  -  A)  may  be  calculated  from  (2.5. 3-2, 3, 4).) 
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Expressions  for  the  coefficients  in  (22)  are  given  below  in  terms  of  the  following  quantities: 

1.  The  equinoctial  elements  ( a,h,k,p,q )  and  the  retrograde  factor  /  (equation  (2. 1.2-2)). 

2.  The  direction  cosines  (q,/3,  7)  of  the  perturbation  symmetry  axis  in  the  equinoctial 
reference  frame  (equations  (2.1. 9-1)). 

3.  The  Kepler  mean  motion  n  (equation  (2. 1.4-3)). 

4.  The  auxiliary  parameter  x  (equation  (2.7. 3-4)). 

5.  The  mean  disturbing  function  U  of  the  perturbation  (equation  (3.1-3)). 

6.  The  coefficients  C1  and  S°  of  the  modified  Fourier  series  expansion  in  L  for  the  short- 
periodic  generating  function  S  (equations  ( 13)— (20)). 

7.  The  cross-derivative  operator  (equation  (2.2-8)). 

8.  The  inclusion  operator  (equation  (2.5.3-18)). 

The  constant  terms  in  the  expansions  (22)  are  given  by: 

2N+1 

C?  =  -  E  {C\p3  +  Si<Tj)  (23) 

j=i 
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The  short-periodic  coefficients  for  the  semimajor  axis  a  are  given  by: 


Cj  =  ~  kkC' + 

+I?~ZU)  [^(i  +  2  )(-hkC’+1  +  ^-=-^S'+2)] 

+i?'‘-J0)[^(i  +  i)(-KV+i  +  tS'+1)] 

+i»-(j)[3^_LXjS)] 

+ir(j)[^u  -  D(^-' + ts'->)j 
+23A,+I(;)[^0  -  2)(MC'-J  +  ^V2)] 

Si  =  I}U)[  ;^(4At/ +  ^-^C' +  MS1)] 

+^(j)|^2*wj 

-^"-3(i)[^0'  +  2)(^^ci+!  +  MS1+2)] 
-^"■20)[^(i  +  i)(‘C'+‘  +  A^+1)] 
-JT'-,y)l^aic>] 

-2)(^-^C2-2  +  MS'-2)] 

A  =  0 
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The  short-periodic  coefficients  for  elements  h  are  given  by: 


si 


-I.!N-2(i)[^0  +  1)C,tI] 

+ir-'u)[^(#.u -vi, ) + + ssM*l 

+I?N(i)[^;(i  -  1)Ci'1j 

ifo)yb(8t,-*c,+fcs')] 

+JJW)[5fcwJ 

+I?Mj)(4^o+2)(*c'+2-‘s,+!)] 

-'**,)  +  jJhfr  - 

+I|JV+,0)[3JS50  - 2)(AC^-2  + 

1  9U  ,  kX  ,  J,  T  ,r  \ 

XnVS^nV^1  /9^7) 


(25) 
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Y-  ! 


:J^J] 


The  short-periodic  coefficients  for  element  k  are  given  by: 

+sw>[sbH,l 

+  2)(-*?«  + 

-If"  *0)[^j(pC.i-,  )  +  xna„5  a/,  -  2nV 

53  =  I-1W[5fe(K‘  +  '‘s')) 

-xr-2o)[^o+i)c3+i] 

-'*4. )  +  + 

+23W+1°)[4^0  -  2)(-Kr’'2  +  **'’)! 

1  dU  hX  , 

'3  = - T~2~aT 2-2\PUiort~lqU,fh) 

\nla 1  ah  n*a* 

-periodic  coefficients  for  element  p  are  given  by: 

ci  -  1 r-u)(1  [ -ci,  +p(c»  -cu -jgj)] 

si  =  ir-'ti) (1  +2P^— -  +p{s -i*  +>C>)1 

D  =  J±±£±£M±l 

Da  2n2a2 

The  short-periodic  coefficients  for  element  9  are  given  by: 

a  =  i?"-o) (l  +££,-■ -  [-/<&  +?(c,i>  -c,i„  -j^)] 

M  -L  r,2  1  nl\vTTI 


d3 

The  short- 


3 

£5  = 


2  n2a2 

(1  +  p2  -f  q2)xIU,cn 
2n2a2 


(26) 


(27) 


(28) 
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The  short-periodic  coefficients  for  the  mean  longitude  A  are  given  by: 


ci  =  -i;o)[ 


X2 


2n2a2(l  +  \) 


(4 hU  -f  — - — C1  +  hkS1)} 


■hkU] 


X2(j  +  2)  /k2-h*. 


SI  = 


(29) 


D«  = 


+I^iSi0x)(hkC,+2 -  ^5,+2)l 

+I,2W  a  +  nvo  +  x) tst  +  4 it)  +  ^(p5,“’  ~,<iS ’’>-■)  -  riws’] 

+^+,°>[2^f^)(AH:,'2+ilT^5"2)] 

2  dU  ,  1  ,9(7  ,  ,01K  X  (  tt  T  Tr  \ 

rfa  da  +  n2a2(l  +  x)  ^  +  *  dk  )  +  n2a2  Vpf/’“7  ) 
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Note  that  the  D,  coefficients  are  simply  related  to  the  first-order  mean  element  rates  Ata 
(given  by  the  right  sides  of  (3.1-1)): 

D,  =  —  (30) 

n 

From  the  central-body  gravitational  potential,  there  are  three  possible  sources  of  explicit 
time-dependence: 

1.  Motion  of  the  central-body  symmetry  axis.  For  the  Earth,  this  is  a  combination  of 
precession  of  the  equinoxes,  nutation,  and  polar  motion. 

2.  Variations  in  the  central-body  rotation  rate. 

3.  Tidal  potential. 

The  principle  effects  of  1  and  2  are  accounted  for  in  SST  by  using  at  each  time  step  the  epoch 
triad  (xg,yB,ZB)  to  evaluate  the  direction  cosines  (u.,/9,-y)  from  (2. 1.9-1)  and  the  rotation 
angle  6  from  (2.7. 1-3,4).  However,  for  the  Earth  the  above  sources  are  thought  to  be  too 
small  or  two  slowly  varying  to  cause  significant  explicit  time-dependence  effects. 

In  order  to  update  the  orbital  elements  in  a  differential  corrections  procedure,  it  is  nec¬ 
essary  to  compute  the  partial  derivatives  with  respect  to  the  mean  elements  of  the  short- 
periodic  variations  (the  ^  in  (2.6-3,  4)).  The  partial  derivatives  of  the  J2  contribution  to 
(22)  are  currently  available  in  the  SST  code,  for  the  special  case  of  zero  eccentricity  and 
replacement  of  (a,/3, 7)  with  the  explicit  formulas  (2. 1.9-3)  in  p  and  q  (thus  motion  of  the 
central- body  symmetry  axis  is  neglected). 


4.2  Third-Body  Gravitational  Potential 

For  a  third-body  point  mass,  the  appropriate  disturbing  function  7£  is  (2.8. 1-6).  From  the 
results  in  Sections  (2.5.1)  or  (2.5.5),  we  can  construct  a  Fourier  series  expansion  in  the  mean 
longitude  A  for  the  first -^rder  short-periodic  variations  r)ia,  as  was  done  by  Green  [1979]. 

However,  for  the  third-body  disturbing  function,  it  is  possible  to  construct  &  finite  Fourier 
series  in  the  eccentric  longitude  F  for  the  qio.  Since  the  D’Alembert  characteristics  are 
bounded,  this  solution  is  of  closed  form  in  the  eccentricity.  In  this  Section  we  outline 
the  construction  of  this  most  desirable  expansion  in  F.  Further  details  can  be  found  in 
[McClain,  1978],  [Cefola  and  McClain,  1978],  [Slutsky  and  McClain,  1981],  and  [Slutsky, 
1983], 

The  disturbing  function  less  its  mean  can  be  written  as 


R  ~  u  =  Re  { it  £ §(2  ~ 6oi)  (wj"  v'~94c.(<*.«  - 


(i) 


j/0 


Here 
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Vn5  =  coefficients  calculated  from  (2.8.2- 1,  2) 

(2„j(7)  =  polynomials  calculated  from  (2. 8. 3-2,  3) 

Cj(a,/?),  Ss(a,/3)  =  polynomials  calculated  from  (2. 5. 3-6) 

Yfs  =  coefficients  of  the  expansion  (2. 7. 1-8) 

The  short-periodic  generating  function  (2. 5. 5-4,  5)  is  easily  obtained  by  integrating  (1): 

{..  N  n  /  n  \n  00  Yr'3e>jX ) 

££E(2 Vn.Qn.lC.(a,0)-iS.(c,f))l  £  -^7— 

K3n=2s=0  j=-oo  lJ  ) 

3*0 


(2) 


The  infinite  series  in  the  mean  longitude  A  in  (2)  can  be  replaced  by  a  finite  series  in  the 
eccentric  longitude  F.  To  see  this,  first  integrate  both  sides  of  (2. 7.1-8)  to  obtain 


oo  ynaeij\  rx  /r\n 

E  =  /  (-)  -  \Yq 

V  J  \aJ 


J  —  —  OO 
3*  0 


(3) 


Next  perform  the  integral  in  (3)  by  using  the  expansion 

0''“=  t  wr*,,F 

vu/  j=—n 


where 


W™  =  Yo~1'’ 


and  by  using  the  change  of  variable  (from  2.5.2-6)) 


The  infinite  series  then  becomes 


j——oo 

3*0 


tj 


dX  =  1 

iddF 

n+1 

Wn+l,s  ijF 

E 

3 

ij 

J=-(n+l) 

3*0 

(4) 


(5) 


(6) 


A) 


(7) 


Replacing  the  infinite  series  in  (2)  with  (7),  recalling  the  relationships  (2. 5. 3-5),  (2.7.1-11), 
(3.1-4),  and  introducing  the  mean  disturbing  function  U  through  (3.2-1),  we  obtain 


5  -  U{F  -  A) 
4- Re 


N  n  /  n  \n  n+* 

£  EB2-M  f)  E 

^3  n=2  ,=0  \tiz*  j=_(n+l) 


[C.(o,S)-tS, 


3*0 


(8) 
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(9) 


The  coefficients  W™  may  be  expressed  in  the  form 

=  e-U-l„,“[Cb-.|(lfc,  h)  -  >sgn(j  -  »)] 

The  functions  w™  possess  the  Jacobi  polynomial  representation 

/l  -2 \ n  (n  +_f)K»  —  -s)! _2\— |a| p|j-d.|j+«l/  \  for  Isl  >  l?l 

«,r=(r^)  (x)  f  |s|£|j|  (10) 

v  +  c '  l(i-  «?)-w/*jrjl-u+*,(x)  for  w  <  ui 

Here  again  x  is  defined  by  (2. 7. 3-4),  e  is  the  orbital  eccentricity,  and 

,  _  e  _  y/h?  +  k2  n  n 

C  l  +  y/T^?  l  +  Vl  —  h2  —  k2  (  ’ 

After  substituting  (9)  into  (8),  we  can  easily  cast  S  into  its  real  form.  The  indices  may 
be  rearranged  as  usual,  and  Kepler’s  equation  (2.1.4)  may  be  used  to  replace  ( F  —  A)  in  (8) 

by 

F  —  A  =  fcsin  F  —  ft  cos  F  (12) 

The  short-periodic  generating  function  can  then  finally  be  written  as  the  finite  Fourier  series 

N+ 1 

S  =  C°  +  U(k  sin  F  —  ft  cosF)  +  (£J  cos  jF  -f  S3  sin  jF)  (13) 

j=i 

Here  the  coefficient  C°  is  (implied  by  <  S  >=  0  and  (2. 5.2-8)) 

C°  =  \c' +h-S'  (14) 

The  coefficients  C 3  and  S3  are 


N  N 


C3  = 


E  E  -f-  {~e-b-*l«,;+1^[sgn(i  -  3)C.(a,/J)V.!iM)  +  Ss(aJ)C\j.3\(k,h)] 

*=0  n=ranx(2,j— 1,»)  7 

+e-U+''«,”+,'*[-C.(o,«Si+.(J:,fc)  +  S,(a,0)C, >,(*,/>)]}  (15) 


N 

S’  =  E 

»=0  n=l 


where 


E  ~{c  b  i|tn”+1’*(C,(a,  P)C\j-t\{k,  ft)  -  sgn (>  -  s)S,(a,  /9)5i,_4|(A;,  A)] 
+e-0+a)u;n+1,i[Ci(Q,  fc)  +  5i(a,  /?)S,+,(fc,  ft)]}  (16) 

Gn.  =  ^(2-^)(^)Vn,(?n4  (17) 
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The  first-order  short-periodic  variations  17 ia  generated  by  the  function  S  given  by  (13)  can 
be  derived  using  equations  (2. 1.8-3),  (2.5.5-10),  and  the  following  (obtained  by  differentiating 
(2. 1.4-2)): 


dF 

dh 

dF 

dk 

dF 

d\ 


(18) 


where  r  is  given  by  (2. 1.4-6).  The  partial  derivatives  of  S  with  respect  to  the  elements 
a,a,/3, 7  follow  by  straightforward  differentiation  of  (13): 


dS 


dc° 


+ 


dU 


d(a,a,0, 7)  d(a,a,{3. 7)  d(a,a,0, 7) 


(k  sin  F  —  h  cos  F) 


N+l 


j=i 


dC> 


,d(a,a,/?,7) 


cos  jF  + 


dSj 


d(a,  a,/?, 7) 


sin  j  F 


(19) 


The  coefficients  and  the  mean  disturbing  function  U  are  functions  of  the  semimajor 

axis  a  through  the  powers  of  the  parallax  factor  alone,  functions  of  the  direction  cosines 
a  and  j3  through  the  polynomials  Ca(a,/3),  St(a,0),  and  functions  of  the  direction  cosine 
7  through  the  polynomials  Qna( 7).  A  finite  Fourier  series  representation  for  §f  may  be 
obtained  by  partial  differentiation  of  (2)  with  the  infinite  series  replaced  by  (2. 7. 1-8),  followed 
by  substitution  of  (4)  and  (9)  and  the  usual  reduction  to  the  real  form 

39  N  . 

+  cos  JF  +  s!x  sin  JF)  (20) 


Here  the  coefficients  are 


r° 


k  .  h 

§Ci  +  2S 


,A 


E  E  G».  h)  -  sgn U  -  *)4(o, «SU-.|(t,  *)I 

*—0  n=max(2,j,a) 


+e-«+*>®-.[C.(o,J8)Ci+.(t.*)  +  S.(o,«Sj+.(t,4)]}  (21) 

Si  =  E  E  G„.{c-U-|)m“[SgnO'-!1)(7.(a,«5|i..|(t,A)  +  5,(Q,(9)Cu..|(t,A)] 

*=0  n=mix(2,j,i) 


+e-«+'V*.[C.(a,i9)5i+.(fc,fc)  -  S.(a,  0)Ci+,(M)]}  (22) 


81 


The  partial  derivatives  of  S  with  respect  to  the  elements  h  and  k  may  be  obtained  by 
differentiation  of  (13)  and  use  of  (18)  and  (20): 


dS  dC°  ,dU  .  _ 
m  =  ~dh+kdKsmF- 


d(hU) 

dh 


4-  cos  F 


(23) 


N+l 

+  E 

j=i 


\%-W-  &*■ )  ™>F + (f  -  k'  -  \s") 


as 


+ 


AT+1 

+  E 

j=i 


a(fct/) 

dk 


+  C°A 


.  „  tdU 
sin  F  —  h-r —  cos  F 
ok 


(24) 


’(?F  -  5s'"  +  5S"‘)  ™’F  +  (f  +  -  5C'") si"'F 


In  the  summations  in  (23)-(24)  Cjx  and  Sjx  are  defined  to  be  zero  for  values  of  j  outside  the 
range  1  <  j  <  N.  The  dependence  of  C\  S>  and  U  on  the  elements  h  and  k  is  through 
the  polynomials  Ct(k,h),S((k,h),  the  eccentricity  e,  and  the  coefficients  w™,Kq*.  In  the 
absence  of  explicit  time-dependence,  the  first-order  short-periodic  variations  can  thereby  be 
written  as  the  finite  Fourier  series 


N+i 

Via  =  c?  +  E  (C?  «»  +  S{  Sin  jF)  (25) 

j=i 

where  the  coefficients  Cf  are  given  by  (2.5.2-15a). 

Green  [1979]  studied  the  effects  of  explicit  time-dependence  by  including  several  terms 
in  the  formulas  (2.5. 1-15)  for  the  coefficients  Cj  and  Sj  of  the  A-expansions  (2. 5. 1-14)  of 
the  short-periodic  variations.  He  used  finite  difference  formulas  to  compute  the  partial 
derivatives  with  respect  to  time  of  the  coefficients  Cj  and  Sj  in  the  A-expansions  (2.5.1-11) 
of  the  osculating  rate  functions.  For  medium  or  high  altitude  Earth  satellites,  he  found 
that  the  Lunar  point  mass  perturbation  varies  quickly  enough  for  explicit  time-dependence 
effects  to  be  significant.  McClain  and  Slutsky  [1988]  also  found  that  the  inclusion  of  explicit 
time-dependence  effects  due  to  the  Moon  and  Sun  improved  the  performance  of  SST  for  high 
altitude  Earth  satellites. 

We  can  include  the  effects  of  explicit  time-dependence  in  the  F-expansions  of  the  present 
Section  by  using  (2.5.1-10)  or  (2.5.5-11)  and  expressions  in  Section  2.5.2.  Specificallly,  the 
first-order  short-periodic  variations  i]ja  including  the  kth  order  time  derivatives  are 

N+k+l 

=  <?  *+  £  (Cl*  cosjF +S!*  sin  jF)  (26) 

J  =  1 


82 


where 


The  recursions  (27)  are  initialized  by 


C°’°  =  Cf 
Cifi  =  C\ 
S{»  =  5/ 


(27) 


(28) 


where  the  C®,  C- ,  Si  are  the  coefficients  in  the  expansion  (25).  The  quantities  U3k,  V*  in  (27) 
are  given  by  the  relations  (2.5.2-12)  with  J  ~  N  +  1  and  m  —  k,  and  X*(j)  is  the  inclusion 
operator  (2.5.3-18). 


4.3  Central-Body  Gravitational  Tesserals 

For  the  central- body  gravitational  tesseral  harmonics,  the  appropriate  disturbing  function 
1Z  is  (2.7.1-16)  with  m  ^  0.  From  the  results  in  Section  2.5.4,  we  can  construct  the  first- 
order  short-periodic  variations  7/,a.  Further  details  beyond  those  given  here  may  be  found  in 
[Proulx,  McClain,  Early,  and  Cefola,  1981]  and  [Proulx,  1981]. 

In  the  absence  of  explicit  time-dependence,  the  first-order  short-periodic  variations  are 
given  by  (2. 5.4-4): 

oo  oo 

Via  =  E  E  lCim  cos(iA  -  m0)  +  sr  sin(jA  -  m6)} 

j=-oo  m=l  V  *  J 


The  C{m  and  Sjm  are  given  by  (2. 5.4-5)  in  terms  of  Fourier  coefficients  C{m  and  Sjm  of 
the  osculating  rate  functions  (2.5. 4-1).  These  latter  coefficients  are  easily  constructed  by 
substituting  the  disturbing  function  (2.7.1-16)  with  m  ^  0  into  (2.2-10).  The  needed  partial 
derivatives  of  1Z  have  the  same  form  as  (3.3-4).  To  illustrate,  we  give  below  the  final  formulas 
for  the  coefficients  C[m  and  S[m: 


ImVZYZKJn-^Pr{GltSnm  - 


"  ImVZ^ZK-n-UPr(GLCnm 


HLCnm) 

+  HLSnm) 


(2) 


(Remember  that  A  is  defined  by  (2.1.6-la).) 

This  theory  has  been  implemented  up  to  a  50x50  gravity  field  model  [Fonte,  1993]. 
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4.4  Atmospheric  Drag 

For  atmospheric  drag,  the  osculating  rate  functions  Fta  are  given  by  (2.2-5)  with  perturbing 
acceleration  q  given  by  (3.4-3).  From  the  results  in  Section  2.5.1,  we  can  construct  an 
expansion  in  the  mean  longitude  A  for  the  first-order  short- periodic  variations: 


Via  =  cos  jA  +  Si  sin; A) 

j=i 


(1) 


where  C-  and  5/  are  given  by  (2.5.1-15)  in  terms  of  Fourier  coefficients  C-  and  S ■  of  the 
A-expansion  (2.5.1-11)  of  the  osculating  rate  functions: 


c;  =  1C  (frq)cosjA<'A 

s‘  = 


(2) 


Of  course,  the  limits  (Ai,A2)  indicate  the  values  of  A  at  atmosphere  entry  and  exit,  and  we 
convert  the  integrals  over  A  into  integrals  over  F  or  L  using  (2.5.2-6)  or  (2.5.3-14).  Green 
[1979]  used  the  A-expansion  (1)  for  atmospheric  drag,  but  indicated  the  desirability  of  an 
expansion  in  another  geometric  variable. 

From  the  results  in  Section  2.5.2,  we  can  construct  an  alternate  expansion  in  the  eccentric 
longitude  F: 

OO 

Via  =  Cf  +  53(C/  cos  jF  +  S\  sin  jF)  (3) 

j= i 

where  C*  and  5/  are  given  by  (2.5.2-14)  in  terms  of  Fourier  coefficients  C,-  and  5/  of  the 
F-expansion  (2.5.2- 1)  of  the  osculating  rate  functions: 


c;  -  lC(^^)cm]FdF 

s:  =  lC(wq)smlF,‘F 


(4) 


But,  as  was  indicated  in  Section  4.1,  the  most  desirable  expansion  for  atmospheric  drag 
is  in  terms  of  the  true  longitude  L.  From  (2.5.3-23) 


K+ 2 


Via  =  Cf  +  £  Am(£  -  A)m  +  cosiL  +  Si  sin  jL) 

J=1 


(5) 


m=l 


where  are  given  by  (2.5.3-24)  (with  the  index  M  =  0)  in  terms  of  the  coefficients 

defined  by 


c/  =  iC{^q)cosjLdL 
sl  =  lC{wq)sinjLdL 


(6) 


D,m  =  0 

From  atmospheric  drag,  there  are  five  possible  sources  of  explicit  time-dependence: 
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1.  Variations  in  the  solar  extreme  ultraviolet  (EUV)  flux. 

2.  Geomagnetic  storms. 

3.  Seasonal  latitudinal  variations  in  the  atmosphere  density. 

4.  Motion  of  the  diurnal  bulge  in  the  atmosphere. 

5.  Motion  of  atmospheric  tides  raised  by  the  Sun  and  Moon. 

Variations  in  the  solar  EUV  flux  can  cause  the  atmospheric  density  at  a  given  altitude  to  vary 
by  up  to  three  orders  of  magnitude  and  are  the  dominant  cause  of  error  in  lifetime  predictions 
for  near-Earth  satellites.  The  primary  source  of  variation  in  the  solar  EUV  flux  is  the  11-year 
sunspot  cycle,  and  it  might  be  assumed  that  variations  with  a  period  of  1 1  years  would  be 
too  slow  to  cause  appreciable  explicit  time- dependence  effects.  The  sunspot  cycle  is  far  from 
sinusoidal,  however.  At  the  beginning  of  each  cycle,  there  is  a  steep  rise  in  solar  EUV  flux 
which  could  conceivably  be  fast  enough  to  cause  significant  explicit  time-dependence  effects. 
In  addition,  there  is  a  secondary  variation  in  solar  EUV  flux  with  a  period  of  about  27  days 
caused  by  the  rotation  of  the  Sun,  which  brings  sunspots  alternately  into  and  out  of  view 
around  the  limb  of  the  Sun.  A  period  of  27  days  is  comparable  to  the  27.3-day  period  of  the 
Lunar  point-mass  perturbation,  which  is  known  to  have  significant  explicit  time-dependence 
effects.  Geomagnetic  storms  can  cause  large  variations  in  the  atmosphere  density  over  time 
scales  of  a  day  or  less.  The  potential  for  significant  explicit  time-dependence  effects  is  obvious. 
The  remaining  sources  of  explicit  time-dependence  effects  for  atmospheric  drag  are  thought 
to  be  too  small  or  too  slowly- varying  to  be  significant.  However,  explicit  time-dependence 
effects  from  atmospheric  drag  have  not  yet  been  studied  with  SST. 

4.5  Solar  Radiation  Pressure 

The  first-order  short-periodic  variations  t/,q  due  to  solar  radiation  pressure  are  formally  iden¬ 
tical  to  the  equations  in  Section  4.4  for  atmospheric  drag,  where  the  perturbing  acceleration 
q  is  given  by  (3.5-6).  For  solar  radiation  pressure,  the  simpler  expansions  (4.4-1)  in  the  mean 
longitude  A  or  (4.4-3)  in  the  eccentric  longitude  F  should  be  adequate.  Also,  Green  [1979] 
found  that  the  explicit  time-dependence  effects  from  solar  radiation  pressure  were  minor. 

As  we  have  seen  in  Section  3.5,  if  the  satellite  is  always  in  sunlight,  the  perturbing 
acceleration  q  can  be  derived  from  a  disturbing  function  K  which  is  nearly  identical  to  the 
third-body  disturbing  function.  Hence  we  can  immediately  obtain  a  finite  Fourier  series  in 
the  eccentric  longitude  F  for  the  rj,a  by  analogy  with  the  results  in  Section  4.2. 


5  Higher- Order  Terms 

Generally,  the  algebraic  complexity  greatly  increases  when  we  try  to  compute  higher-order 
mean  element  rates  and  short-periodic  variations.  Also,  it  is  assumed  that  higher-order 
terms  due  to  most  perturbations  are  usually  negligible.  In  this  chapter  we  report  on  those 
few  higher-order  effects  which  have  been  studied  to  date  with  SST. 
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5.1  Second-Order  Amp  and  r)iap  Due  to  Gravitational  Zonals  and 
Atmospheric  Drag 

The  second-order  mean  element  rates  Aial 3  and  short-period  variations  qlQp  due  to  two  per¬ 
turbations  expanded  in  A  may  be  constructed  as  shown  in  Section  2.5.6.  We  can  calculate 
analytically  the  Fourier  coefficients  and  Sfj  of  the  expansions  in  A  for  the  osculating 
rate  functions  Fn  due  to  the  central-body  gravitational  zonal  harmonics  by  substituting  the 
disturbing  function  (2.7.1-16)  with  m  =  0  into  equations  (2.2-10).  We  can  calculate  nu¬ 
merically  the  Fourier  coefficients  C/2  and  S-2  of  the  expansions  in  A  for  the  osculating  rate 

dS* 

functions  Fn  due  to  atmospheric  drag  from  (4.4-2).  The  partial  derivatives  and 
needed  in  (2.5. 6-5)  can  be  calculated  analytically  for  the  central-body  gravitational  zonals 
and  by  numerical  quadrature  for  atmospheric  drag. 

At  the  present  time,  the  only  terms  in  these  analytical  formulas  which  are  available  in 
the  SST  code  are: 

1.  The  J^-squared  auto-coupling  mean  element  rates  Am,  correct  to  first  power  of  the 
eccentricity  and  with  (0,^,7)  replaced  by  the  explicit  formulas  (2. 1.9-3)  in  p  and  q. 
These  terms  were  constructed  with  the  MACSYMA  algebraic  utilities  of  [Zeis,  1978] 
and  [Bobick,  1981]. 

2.  The  J2-  squared  auto-coupling  short-periodic  variations  r/,u,  correct  to  zero  power  of 
the  eccentricity  and  with  (o,/3, 7)  replaced  by  the  explicit  formulas  (2. 1.9-3)  in  p  and 
q.  These  terms  were  constructed  in  [Zeis,  1978]  and  corrected  in  [Green,  1979]. 

3.  The  J2I drag  cross-coupling  mean  elements  rates  Am,  correct  to  zero  power  of  the 
eccentricity  and  with  the  retrograde  factor  I  =  1  and  (a,/?, 7)  replaced  by  the  explicit 
formulas  (2. 1.9-3)  in  p  and  q.  These  terms  were  constructed  in  [Green,  1979]. 

Green  [1979]  studied  the  second-order  effects  of  J2  and  drag  using  a  combination  of 
analytical  and  numerical  methods.  He  found  that  the  drag/  J2  cross-coupling  terms  Ai2i  cause 
significant  effects  for  low  altitude  satellites.  He  found  that  Izsak’s  J2  height  correction  (1) 
applied  to  the  density  determination  in  the  formulas  (3.4-1,  3,  4)  gave  a  good  approximation 
to  the  An  1  terms.  The  following  expression  added  to  the  height  H  is  the  J2  short-periodic 
correction  (from  Section  4.1)  to  the  radial  distance  r: 


A  r  = 


J2R 2 

4(1  —  e2)a 


sin2  i  cos  2 (/  +  w)  -(-  (3  sin2  i  —  2) 


ecos  / 

1  +  y/l  ~  e2  + 


2yT=~e2  \ 
1  +  e  cos  /  / 


(1) 


Green’s  explanation  of  why  Izsak’s  J2  height  correction  works  is  that  adding  the  J2  short 
periodics  to  the  elements  in  the  drag  osculating  rate  functions  and  then  averaging  is  equiva¬ 
lent  to  adding  the  drag/  J2  cross-coupling  terms  to  the  first-order  mean  element  rates,  if  we 
neglect  products  of  short-periodic  variations: 


<  F i2(®i  +  Vni  +  »?6 i)  > 


d  F 

~  <  Fa(au-  ■•a6)  >  +X!  <  0^(01  ,--’ae)Vji  > 
~  Ai2  +  Ai  21 


(2) 
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5.2  Second-Order  r ]iQfj  Cross- Coupling  Between  Secular  Gravita¬ 
tional  Zonals  and  Tesseral  Harmonics 


In  high-order  shallow  resonance  orbits,  the  tesseral  harmonics  which  contribute  the  most 
significant  short-period  motion  are  likely  to  be  those  with  degree  and  order  centered  around 
the  resonant  order.  For  such  orbits,  the  second-order  short- periodic  variations  due  to  cross¬ 
coupling  between  these  tesseral  harmonics  and  the  J2  secular  terms  may  also  be  significant. 
In  this  Section  we  outline  how  to  construct  these  critical  short  periodics.  For  further  details 
and  a  discussion  of  numerical  results  see  [Cefola,  1981]  and  [Cefola  and  Proulx,  1991]. 

For  the  present  purpose,  we  retain  in  the  expansions  (2.5.1-10)  of  the  osculating  rate 
functions  Fn  due  to  the  central- body  gravitational  zonal  harmonics  only  the  mean  element 
rates  An  given  by  (3.1-1): 

Ftl  %  Ati  (1) 

Furthermore,  we  completely  neglect  the  first-order  short-periodic  variations  r/tl  due  to  the 
zonal  harmonics: 

Vn  ~  0  (2) 

The  osculating  rate  functions  Fa  due  to  the  tesseral  harmonics  may  be  expanded  in  the 
Fourier  series 

Fn  =  T}dm  cos(jA  —  mO)  +  Sjm  sin(jA  —  m&)\  (3) 

j,m 

The  first-order  short-periodic  variations  due  to  the  tesseral  harmonics  are  then  given  by 
(2. 5.4-4): 

fh  2  =  ”  COS0'A  ~  m&)  +  Si™  sin  O' A  -  mO)}  (4) 

where,  in  the  absence  of  explicit  time-dependence,  the  Cf m  and  SJ,m  are  related  to  the  Fourier 
coefficients  C-m  and  Sjm  by  (2.5. 4-5). 

To  obta:;i  the  second-order  cross-coupling  terms,  we  need  to  construct  the  functions 
G{ i2  +  G{2i  from  (2.3-27): 


G,  12  +  Gi2i 


15n 

+  -r-^BieVuVn 


Substituting  (l)-(4)  into  (5)  yields 

(7,12  +  Gm  «  ^2  [cr  cosO’A  —  mO)  -f  3; m  sinO’A  —  mO)] 


where 


cr 


sr 


-jA61Sr 

+  jAeiC?m 


(5) 


(6) 


(7) 
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r 


The  cross-coupling  short- periodics  are  then 

+  i?, 2i  =  J2  [C-mcos(jA  -  m9)  +  sin(;A  -  m0)]  (8) 

where  the  coefficients  C3m  and  are  given  in  terms  of  C\  and  S3t  by  the  relations 
(2. 5.4-5).  The  partials  of  the  J2  mean  element  rates  needed  in  (7)  are  given  by  equations 

(3.1-12).  The  partials  and  -g£—  of  the  tesseral  first-order  short-periodic  coefficients  are 

related  as  usual  to  the  partials  and  of  the  tesseral  osculating  ra'e  coefficients, 

which  may  be  obtained  by  differentiating  formulas  such  as  (4.3-2). 

Code  based  on  this  approximate  theory  has  been  developed  only  for  the  J2  secular /m- 
daily  terms  (j  =  0  in  (4)),  with  finite  differences  used  to  obtain  the  partials  of  the  m-daily 
coefficients.  Construction  and  programming  of  a  complete  second-order  theory  for  a  double- 
averaged  perturbation  expanded  in  A,0  has  yet  to  be  accomplished. 


6  Numerical  Methods 

The  numerical  methods  which  are  currently  used  in  SST  are  standard.  In  this  chapter  we 
record  the  essential  mathematical  formulas.  Further  details  may  be  found  in  any  numerical 
analysis  textbook  (e.  g.,  [Ferziger,  1981)). 


6.1  Numerical  Solution  of  Kepler’s  Equation 

The  equinoctial  form  of  Kepler’s  Equation  is  (2. 1.4-2): 

A  =  F  +  h  cos  F  —  k  sin  F 

This  equation  can  be  solved  iteratively  using  Newton’s  method: 

F0  =  A 

^  r,  (Ft  +  h  cos  F,  -  k  sin  -  X\  . 

Fi+ 1  =  Fi  -  — - = — 7 - p —  for  *  =  0, 1,2, 

\  1  —  h  sin  Fi  —  k  cos  Fi  ) 


(1) 

(2) 


6.2  Numerical  Differentiation 

We  need  to  differentiate  functions  in  order  to  obtain  the  mean  element  rates,  short-periodic 
variations,  and  partial  derivatives  for  state  estimation.  Analytical  formulas  are  preferable  if 
possible  to  obtain,  because  of  their  greater  precision.  However,  the  derivatives  of  a  function 
can  be  approximated  by  finite  difference  schemes. 

We  suppose  /( i)  is  a  smooth  function  of  x.  Then  the  central  difference  approximation 
for  the  derivative  of  f(x)  is 

a, 
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The  error  in  this  approximation  is 


df,  , 

rJx)- 


f{x  +  A)  -  f{x  -  A) 

2A 


-  dPf 

6  dx 3 


(£),  i  —  A  <  £  <  a:  +  A 


For  example,  Green  [1979]  used  the  central  difference  approximation  (1)  to  calculate  the 
partial  derivatives  needed  for  state  estimation  (see  Section  2.6).  He  obtained  good  results 
with  a  step  size  of  A  =  10-5x,  using  double  precision. 


6.3  Numerical  Quadrature 

We  need  to  integrate  functions  in  order  to  obtain  the  mean  element  rates  and  short-periodic 
coefficients.  Analytical  formulas  are  preferable  if  possible  to  obtain,  because  of  their  greater 
precision.  Also,  they  are  more  computationally  efficient  since  analytical  formulas  need  to 
be  evaluated  only  once  per  mean  equations  integration  step,  whereas  numerical  integration 
requires  evaluation  at  each  abscissa  of  the  quadrature  [Long  and  McClain,  1976].  However, 
numerical  evaluation  of  integrals  of  the  type 

/  f{x)dx  (1) 

J  a 

is  mandatory  for  the  computation  of  the  mean  element  rates  and  short-periodic  coefficients 
involving  atmospheric  drag  or  solar  radiation  pressure  with  eclipsing.  Since  the  substitution 

e  2x  -  (a  +  b) 

*  “  b^a 

tranforms  the  integral  ( 1 )  into 

=  ^ /_'/(«)<({  (3) 

we  can  restrict  our  discussion  to  integrals  with  limits  between  —1  to  +1  without  loss  of 
generality. 

A  quadrature  formula  approximates  an  integral  by  a  weighted  sum  of  the  values  of  the 
integrand  at  points  on  the  interval  of  integration: 

f  f(t)dt  »  £  w,m),  -1  <  6  <  6  <  •  •  -6.  <  1  (4) 

.= 1 

An  evaluation  of  different  quadrature  formulas  has  shown  the  Gaussian  quadrature  formulas 
to  be  generally  efficient  [Early,  1975].  The  weight  factors  w,  for  Gaussian  quadratures  have 
been  tabulated,  and  the  abscissas  &  are  simply  the  zeroes  of  the  Legendre  polynomial  of 
degree  n.  The  error  in  the  Gaussian  quadrature  formula  is 


(2) 


>=i 


22n+1(n!)4  <Pnf{  0 
(2n  +  l)(2n!)3  d £2"  ’ 


-]<£<] 


(5) 


A  polynomial  of  degree  2n  —  1  is  integrated  exactly. 

The  appropriate  number  n  of  abscissas  in  the  Gaussian  quadrature  formulas  needed  for 
SST  can  vary  from  12  to  96,  depending  on  the  highest  frequency  components  contained  in 
the  function  to  be  integrated.  For  example,  Green  [1979]  found  that  if  the  first  10  pairs  of 
short-periodic  coefficients  are  to  be  retained  in  (2.5.1-13),  the  number  n  for  the  integrals 
(2.5.1-11)  must  be  at  least  48. 


6.4  Numerical  Integration  of  Mean  Equations 

The  averaged  equations  of  motion  (1-2)  may  be  solved  with  a  Runge-Kutta  numerical  inte¬ 
gration  method.  We  consider  the  following  system  of  ordinary  differential  equations: 

(1) 

Here  x  denotes  the  column  matrix  of  mean  elements,  and  f  denotes  the  column  matrix  of 
mean  element  rates.  We  divide  the  t- axis  into  points  (<i,<2»’“)  of  equal  width  ft,  and  let 
Xi  =  x(tt).  Then  the  standard  fourth-order  Runge  Kutta  algorithm  is 


x«+i  —  x,-  +  —  (Wj  +  21<2  +  2k3  +  k4)  for  i  =  1,2,3,*  •  •  (2) 

o 

where 

kj  =  f(x„t,) 

k2  =  f(x,  +  —  kj,f,  +  — ) 

2  2  (3) 

A  A  K 

k3  =  f(x,  +  yk2,f,  +  -) 

let  =  f(x,  +  Ak3,*i  +  A) 

The  error  in  the  formulas  (2)  is  bounded  by 


CA 


s£xi 

dt 5 


(4) 


where  C  is  a  constant. 

Since  the  mean  element  rates  depend  only  on  slowly  varying  quantities,  step  sizes  A 
of  a  day  or  more  can  usually  be  used.  The  integrator  time  step  A  should  be  |  or  less  of 
the  minimum  period  r  of  the  oscillations  included  in  the  mean  equations  of  motion.  Some 
limitations  are  the  period  of  orbital  precession  due  to  J2  and  the  period  of  the  moon. 

Initial  values  of  the  mean  elements  a,(*i)  can  be  obtained  from  initial  values  of  the 
osculating  elements  dj(<j)  by  either  of  two  methods: 


1.  Numerically  integrate  the  VOP  equations  of  motion  over  a  time  interval  at  least  as 
long  as  the  period  of  the  largest  significant  short-periodic  effect  (usually  one  or  two 
satellite  orbits  -  see  [McClain  and  Slutsky,  1980]),  and  then  use  a  differential  correction 
procedure  to  find  the  initial  mean  elements  which  give  the  best  least-squares  fit  between 
the  SST  trajectory  and  the  Cowell  trajectory. 
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2.  Use  successive  substitution  into  the  near-identity  transformation  (1-1)  until  a  specified 
agreement  is  reached: 

«?(<l)  =  M*  l) 

a?+1(*i)  =  -T]i[akl(t1),---a^{t1),tl]  for  k  =  0,1,2,- •• 

This  method  is  faster  than  method  1 ,  but  may  require  the  inclusion  of  a  comprehensive 
set  of  short-periodic  variations  to  avoid  a  large  bias  in  the  initial  mean  elements. 

It  should  be  pointed  out  that  the  time  averages  of  the  osculating  elements  over  some  time 
interval  are  generally  not  a  good  approximation  to  the  mean  elements  [Early,  1986]. 

6.5  Interpolation 

Since  the  mean  elements  and  short-periodic  coefficients  are  slowly  varying,  their  values  at 
desired  times  not  coinciding  with  the  mean  equation  step  times  can  be  computed  by  relatively 
low  order  interpolation  formulas. 

First,  suppose  that  at  distinct  times  (tj, •  •  •  tn)  we  know  the  values  [/(<i),  •  •  •  f{tn)]  of  a 
smooth  function  f(t).  In  Lagrange  interpolation  we  approximate  f(t)  by  a  polynomial  of 
degree  n  —  1  passing  through  the  known  values: 

/W«E /((.)£.  W 

1=1 

Here 

L  (t)  =  (*~  *0  •••(<-<,-!)(< •••(<-*») 

,U  (U-ti)---  ( t ,  -  tt-r)(ti  -  tm)  •  •  •  (t,  -  tn) 

Note  that 

Li(tj)  =  6ij 

The  error  in  the  Lagrange  interpolation  formula  is 

m - 1, mm)  = (t <.<«<<.  w 

Lagrange  interpolation  is  currently  used  to  interpolate  the  short-periodic  coefficients,  the 
velocity  vector,  and  the  partial  derivatives  needed  for  differential  correction.  An  adequate 
order  n  --  1  has  been  found  to  be  3  (4  interpolator  points)  [Taylor,  1978]. 

Next,  suppose  that  at  distinct  times  (ti,  •  •  •  tn)  we  know  both  the  values  [/(<i),  ■  ■  •  /(<n)] 
and  the  derivatives  [/(<i),  •  •  •  /(^n)]  of  a  smooth  function  f(t).  In  Hermite  interpolation  we 
approximate  f(t)  by  a  polynomial  of  degree  2n  —  1  passing  through  the  known  values  and 
derivatives: 


0) 

(2) 

(3) 


m  «  E  {[i  -  2(t  -  *o^te)][MO]2/te) + (<  -  ti)[Li(t)]2m} 

i= 1 


(5) 
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Here  again  £,,-(<)  axe  the  Lagrange  basis  functions  (2).  The  error  in  the  Hermite  interpolation 
formula  is 


[( t-u)-jt-tn))*dhf 
(2n)!  d( 


U<Z<tn 


(6) 


Hermite  interpolation  is  currently  used  to  interpolate  the  mean  elements  and  the  position 
vector.  An  adequate  order  2n  —  1  has  been  found  to  be  5  (3  interpolator  points). 
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